Image model based on n-pixels and defined in algebraic topology, and applications thereof

ABSTRACT

A computational image model comprises an image support including a structure of n-pixels comprising pixel faces, quantities related to image features, and an algebraic structure relating the quantities to the n-pixels and/or pixel faces, the algebraic structure comprising algebraic operations defining a relation between the quantities. A method of computationally modelling an image comprises producing an image support including a structure of n-pixels comprising pixel faces, defining quantities related to image features, and relating the quantities to the n-pixels and/or pixel faces through an algebraic structure, and relating the quantities to each other through algebraic operations.

FIELD OF THE INVENTION

The present invention relates to an image model based on n-pixels. Morespecifically, the present invention is concerned with an image modelbased on n-pixels, defined in algebraic form and applicable, inparticular but not exclusively, for the resolution of diffusion andoptical flow, and for the deformation of curves.

BACKGROUND OF THE INVENTION

People have a notion of what an image is. For instance, forpsychologists the image is linked to the shape of objects, their depth,the relationship of these shapes and their perceptual organization.

Artists are focused on how features such as shape, color, andperspective are organized to represent a scene that may originate intheir imagination.

Physicists are concerned with the physical phenomena produced by a givenscene and how they are represented in the image.

Neurophysiologists regard images through visual phenomena in humans andanimals, such as contrast sensitivity, Mach bands, contrast, constancy,and depth perception.

While a precise definition of the image is elusive, it is clear that forcertain people an image is the visualization of physiological,perceptual, or physical phenomena and for others it is related tosemantic content. Whatever the term image means, the intention is toestablish a foundation upon which images of all forms and contents canbe discussed with minimal confusion.

In image processing and computer vision, the foundation is related tothe understanding of the image formation process. Light generated by asource is modified by the objects of the scene. The modified light iscaptured by a system acquisition device, transformed into an appropriateform and displayed on a physical support. The content of the resultingimage and, consequently, its further use, both deperid on the propertiesof the light (structured or not, spectral properties such as the rangeof wavelengths, the number of ranges and the intensity), thecharacteristics of the acquisition device, the transformation(analog-to-digital, pre-processing), the display format (temporal orspatial organization of image elements, film, vector, raster). From theautomatic processing point of view, the image is enhanced to improve itsperceptual quality or to make some of its features explicit. Usually,its content is analyzed via a successive reduction process to constructa more descriptive representation in terms of relevant features, whichcan be used more effectively by a decision system (car, robot, etc.) orto help human beings in their daily tasks.

One of the concerns here is to focus on the data structure for imagesand its consequences on the processing scheme.

Algebraic topology concepts are a key to representing images. Algebraictopology is well-known domain in mathematics [10, 6, 9]. The literatureof algebraic topology offers wide knowledge that can be used for images.The specialists for algebraic topology, however, have made no effort toimplement their knowledge into computer vision and image processing.Instead of using algebraic topology, the specialists of image processingand computer vision have limited themselves to develop solutions basedon the topology and on discrete geometry [7, 8].

Algebraic topology was first introduced into image processing andcomputer vision by Allili and Ziou for topological feature extractionand shape representation in binary images [1, 2, 15]. This technology isused by Allili, Mischaikow and Tannenbaum [3] in medical binary images.Auclair et al. [4] also used algebraic topology for linear andnon-linear isotropic diffusion as well as for optical flow in gray levelimages. P. Poulin et al. [12] used algebraic topology for snakes andelastic matching in gray level images.

In image processing and computer vision, several image models have beenaccepted and are in recurrent use since several decades. These imagemodels integrate both the image support and the local quantitiesassociated with this support. The image support is formed by pixels.With each pixel is associated a scalar quantity called a gray level, ora vector quantity called either color at the perceptual level ormultispectral at the signal level. Existing models differ primarily interms of how the image support (definition and the organization ofpixels) and how values are formulated. The well-known image models arethe function, the random process, and the ordered set. The image is afunction L_(x)×L_(y)→G^(m), where L_(x)={1, . . . ,N_(x)} and L_(y)={1,. . . , N_(y)}, N_(x)×N_(y) is the resolution of the image, G={1, . . ., n}, where n is the maximal quantity and m the number of image bands.In the case of a binary image n=2, image processing has roots in graphtheory, language theory, logic, discrete geometry, and so on. If n>1,usually the image is modelled as a real function (analogue image whereG,L_(x),L_(y) ⊂ R). In this case, function theory, functional analysis,catastrophe theory, differential equations, and differential geometryconstitute the foundation. An image can also be modelled as a collectionof random variables X(i,j)|(i,j) ε L_(x)×L_(y). In this case, theprobability density function, moments, sufficient statistics, timeseries, and the Markov processes are the roots. When the image ismodelled as an ordered set, discrete mathematics and mathematicalmorphology are the foundation.

Since the introduction of mathematical morphology, efforts ofresearchers in these fields have been focused on the use of more andmore complex mathematical, physical or computer concepts as formalism ofspecific problems (edge detection, image segmentation, optical flow, anddeformation) without questioning the image model. The definition of anew image model can lead algorithms that are designed on new basis. Animage is a physical or mathematical quantity where variables (imagesupport) represent geometrical or temporal elements such as points,lines, surfaces, and times. For example, the image, as a functionL_(x)×L_(y)→G^(m), can be defined by both the geometrical andtopological properties of the domains L_(x)×L_(y) and the topological,geometrical and analytical properties of G^(m). Although existing imagemodels have deep roots in mathematics or in physics, the variables, thequantity and the association between them are not well-defined. For agiven computer vision or image processing task, no formal mechanism isgiven for the integration of physical, topological, geometricalproperties of objects and their behaviours as a part of the image model.Consequently, the resulting computational schemes are non-modular andsometimes not easy to reproduce.

SUMMARY OF THE INVENTION

According to the present invention, there is provided a computationalimage model, comprising an image support including a structure ofn-pixels comprising pixel faces, quantities related to image features,and an algebraic structure relating the quantities to the n-pixelsand/or pixel faces, the algebraic structure comprising algebraicoperations defining a relation between the quantities.

The present invention also relates to a method of computationallymodelling an image, comprising producing an image support including astructure of n-pixels comprising pixel faces, defining quantitiesrelated to image features, and relating the quantities to the n-pixelsand/or pixel faces through an algebraic structure, and relating thequantities to each other through algebraic operations.

Still in accordance with the present invention, there is provided acomputational framework for solving a heat transfer problem, comprising:

-   -   producing an image support including a structure of n-pixels,        the image support comprising:        -   q-pixels respectively translating the n-pixel algebraically,            wherein q ε {1, 2, . . . , n}, and wherein each q-pixel            includes (q−1)-faces, (q−2)-faces, . . . , (q−q)-faces;        -   geometrical complexes each being a collection of q-pixels;        -   q-chains respectively expressing the geometrical complexes            in algebraic form, each q-chain being a linear combination            of all the q-pixels of the geometrical complex;        -   in the geometrical complexes, q-cochains which are relations            associating quantities related to image features to the            q-pixels and/or faces of said q-pixels; and        -   a coboundary defining a relation between q-cochains;    -   computing a q-cochain T of a first geometrical complex as the        location of unknown temperatures;    -   computing a q-cochain H of the first geometrical complex as a        global temperature variation;    -   finding a q-cochain ε of a second geometrical complex as a        global energy variation, as a function of the q-cochain H        through a linear transformation;    -   finding the q-cochain ε as a function of the q-cochain T;    -   defining a q-cochain G of the first geometrical complex from the        q-cochain T-through a first coboundary operation, transforming        the q-cochain G into a q-cochain Q of the second geometrical        complex, and defining, from the q-cochain Q and through a second        coboundary operation, a q-cochain D of the second geometrical        complex as a global diffusion;    -   defining a q-cochain S of the second geometrical complex as a        global source;    -   establishing a relation between the q-cochains ε, D and S.

The present invention further relates to a computational framework fortwo-dimensional active contour model, comprising:

-   -   producing an image support including a structure of n-pixels,        the image support comprising:        -   q-pixels respectively translating the n-pixel algebraically,            wherein q ε {1, 2, . . . , n}, and wherein each q-pixel            includes (q−1)-faces, (q−2)-faces, . . . , (q−q)-faces;        -   geometrical complexes each being a collection of q-pixels;        -   q-chains respectively expressing the geometrical complexes            in algebraic form, each q-chain being a linear combination            of all the q-pixels of the geometrical complex;        -   in the geometrical complexes, q-cochains which are relations            associating quantities related to image features to the            q-pixels and/or faces of said q-pixels; and        -   a coboundary defining a relation between q-cochains;    -   computing a displacement q-cochain D of a first geometrical        complex;    -   computing a strain q-cochain S of a second geometrical complex,        comprising:        -   defining an approximate strain function {tilde over (ε)}(x)            as a function of the q-cochain D;        -   expressing the q-cochain S as a function of the approximate            strain function and relative positions of the first and            second geometrical complexes; and    -   computing a force q-cochain F of the second geometrical complex        as a coboundary of the strain q-cochain S.

The foregoing and other objects, advantages and features of the presentinvention will become more apparent upon reading of the followingnon-restrictive description of illustrative embodiments thereof, givenby way of example only with reference to the accompanying drawings.

The present specification refers to various references. These referencesare herein incorporated by reference in their entirety.

BRIEF DESCRIPTION OF THE DRAWINGS

In the appended drawings:

FIG. 1 is an example of a 2-cube; the orientation is given by definition[0,1]×[0,1];

FIG. 2 a is an example of subdivision that is a cubical complex, andFIG. 2 b is another example of subdivision that is not a cubicalcomplex;

FIG. 3 illustrates, in solid lines, an example of a primary cubicalcomplex and, in dashed lines, an example of a secondary cubical complex;

FIG. 4 a is an example of original image and FIG. 4 b is an example of aresulting smoothed image;

FIGS. 5 a, 5 b and 5 c illustrate a body in 3D space at time t. In FIG.5 a, x₁ is the location of the particle X₁, n is a vector normal to thesurface, dS is an infinitesimal surface patch and dV is an infinitesimalamount of volume. In FIG. 5 b, f_(e) is an external force applied on dS,ρ is the mass and b is a body force applied on dV. In FIG. 5 c, q is theheat flow passing through dS and r is the heat produced by dV;

FIGS. 6 a, 6 b and 6 c are three examples of q-pixels in R² (I₁×I₂).FIG. 6 a illustrates an example of 0-pixel where I₁={2} and I₂={1}, FIG.6 b an example of 1-pixel where I₁=[2,3] and I₂={1}, and FIG. 6 c anexample of 2-pixel where I₁=[2,3] and I₂[1,2];

FIG. 7 is an example of coboundary operation;

FIG. 8 a is an example of projection onto the tangential part of thedomain, and FIG. 8 b is an example of projection onto the normal part ofthe domain;

FIGS. 9 a and 9 b are examples of cochains for one 3-pixel of K^(p1).FIG. 9 a illustrates an example of cochain T while FIG. 9 b illustratesan example of cochain H;

FIG. 10 a is an example of cochain Q for one 3-pixel of K^(s), and FIG.10 b is an example of cochain D for one 3-pixel of K^(s);

FIG. 11 a is an example of cochain T for one 3-pixel of K^(p), and FIG.11 b is an example of cochain G for one 3-pixel of K^(p);

FIG. 12 is an example of three different paths between two points;

FIG. 13 is an example of computational scheme for an unsteady problemwith no source;

FIG. 14 is an example of two cubical complexes for a 5×5 image;

FIG. 15 is an example of γ_(p) for one 2-pixel of K^(p);

FIG. 16 a is an example of γ_(E) (in dashed lines) for one 2-pixel of K³intersecting four pixels of K^(p), FIG. 16 b is an example of cochains Tand G, and FIG. 16 c is an example of cochain Q;

FIG. 17 is an example of one 3-pixel of K^(s′) surrounding one 1-pixelof K^(p′);

FIG. 18 is an example of λs for one 2-pixel of K^(p);

FIGS. 19 a, 19 b and 19 c are examples of one 3-pixel of K^(s)intersecting four 3-pixels of K^(p), for cochain T (FIG. 19 a), cochainG (FIG. 19 b), and cochain Q (FIG. 19 c);

FIG. 20 is an example of two 2-pixels of K^(s) with λ's;

FIG. 21 a is an example of physics-based isotropic diffusion σ={2.0,4.0, 5.0}, and FIG. 21 b is the example of physics-based isotropicdiffusion of FIG. 21 a convolved, for same scales;

FIGS. 22 a, 22 b and 22 c are examples of first images of threesequences, a rotating sphere sequence (FIG. 22 a) a Hamburg taxisequence (FIG. 22 b) and a tree sequence (FIG. 22 c);

FIG. 23 a is an example of flow pattern computed for the rotating spheresequence of FIG. 22 a using a physics-based method, and FIG. 23 b is anexample of flow pattern computed for the rotating sphere sequence ofFIG. 22 a using an iterative finite difference method;

FIG. 24 a is an example of flow pattern computed for the Hamburg taxisequence of FIG. 22 b using the physics-based method, and FIG. 24 b isan example of flow pattern computed for the Hamburg taxi sequence ofFIG. 22 b using the iterative finite difference method;

FIG. 25 a is an example of flow pattern computed for the tree sequenceof FIG. 22 c using the physics-based method, and FIG. 25 b is an exampleof flow pattern computed for the tree sequence of FIG. 22 c using theiterative finite difference method;

FIG. 26 a is an example of a first image of the tree sequence of FIG. 22c with white noise added (standard deviation of 10);

FIG. 27 a is an example of flow pattern computed for the tree sequenceof FIG. 22 c with white noise added using the physics-based method, andFIG. 27 b is an example of flow pattern computed for the tree sequenceof FIG. 22 c with white noise added using the iterative finitedifference method;

FIG. 28 a is a section of the peppers image (σ=5) of FIG. 21 acorresponding to the original image with noise added, FIG. 28 b is asection of the peppers image (σ=5) of FIG. 21 a obtained with the PBmethod, and FIG. 28 c is a section of the peppers image (σ=5) of FIG. 21a obtained with the FD method;

FIGS. 29 a, 29 b, 29 c and 29 d are examples of a section of the peppersimage (σ=5) of FIG. 21 a corresponding to the original image (FIG. 29a), corresponding to the original with white noise added (FIG. 29 b,obtained using the PB method (FIG. 29 c), and obtained with the FDmethod (FIG. 29 d);

FIG. 30 a is an example of PB method for σ={1.0, 3.0, 5.0}, and FIG. 30b is an example of FD method for the same scales;

FIG. 31 a is an example of original image, and FIG. 31 b is an exampleof image with noise added (standard deviation=10);

FIG. 32 a is an example of PB method for σ={4.0, 8.0}, and FIG. 32 b isan example of FD method for the same scales;

FIG. 33 a is an example of PB method, and FIG. 33 b is an example of FDmethod with σ=8.0;

FIG. 34 is an example of a body of arbitrary size, shape and material,where ΔS is a surface patch, f is a vector of surface forces applied onΔS, ΔV is an amount of volume with a mass ρ, and b is a vector of bodyforces applied on ΔV;

FIG. 35 is an example of cutting plane passing through a point andpartitioning the body into two sections;

FIG. 36 is an example of a force Δf acting on a cutting plane withnormal vector n;

FIG. 37 a is an example of stresses on the original body, FIG. 37 b isan example of normal stress after an extension of the body, FIG. 37 c isan example of normal stress after a compression of the body, and FIG. 37d is an example of shear stress after a distortion of the body;

FIG. 38 is an example of the component of the stress in the direction ofx_(I);

FIGS. 39 a and 39 b are examples of the deformation (FIG. 39 a) anddistortion (FIG. 39 b) of a body subjected to stresses; in both FIGS. 39a and 39 b, the rectangle ABCD has been deformed or sheared into A′B′CD;

FIG. 40 is an example of normal strain of a body;

FIG. 41 is an example of shear strain in a body;

FIG. 42 is an example of a body B in motion and subjected to forces,wherein t_(i) ^(n) and pb_(i) are respectively the traction and bodyforces in the direction of x_(I);

FIG. 43 a is an example of kinematic equation, FIG. 43 b is an exampleof constitutive equation, and FIG. 43 c is an example of conservationequation;

FIG. 44 is an example of decomposition of the linear elasticity probleminto basic laws;

FIGS. 45 a, 45 b and 45 c are three examples of q-pixels in R² (I₁×I₂).More specifically, FIG. 45 a is an example of 0-pixel for I₁={2} andI₂={1}, FIG. 45 b is an example of 1-pixel for I₁=[2,3] and I₂={1}, andFIG. 45 c is an example of 2-pixel for I₁=[2,3] and I₂=[1,2];

FIG. 46 is an example of the coboundary operation;

FIG. 47 a is an example of cochain U for a 3-pixel of K^(p), and FIG. 47b is an example of cochain D for a 3-pixel of K^(p);

FIGS. 48 a and 48 b are example of cochains S and F, respectively, for a3-pixel of K^(p);

FIG. 49 is an example of two cubical complexes for a 5×5 image;

FIG. 50 is an example of a 2-pixel of K^(p) and the topologicalquantities associated with it;

FIG. 51 a is an example of γ_(F) in dashed lines, FIG. 51 b is anexample of 2-cochains D and U, and FIG. 51 c is an example of 2-cochainS;

FIG. 52 is an example of five adjacent vertices in a curve and itsdeformation;

FIG. 53 is a table that shows typical values of the Poisson's ratio andYoung's modulus of elasticity for some materials;

FIG. 54 a is an example of initial curve, and FIG. 54 b is a bright lineplausibility image;

FIG. 55 a is an example of results obtained for an aerial image usingthe PB method, and FIG. 55 b is an example of results obtained for anaerial image using the FEM method;

FIG. 56 is an example of initial curve for a SAR image;

FIG. 57 is an example of line plausibility for a SAR image;

FIG. 58 is an example of road correction for a SAR image with the PBmethod;

FIG. 59 is an example of road correction for a SAR image with the FEMmethod;

FIG. 60 is an example of initial curve;

FIG. 61 is an example of bright line plausibility image;

FIG. 62 is an example of corrections for a multiband image (PB method);

FIG. 63 is an example of corrections for a Landsat 7 image (FEM);

FIG. 64 a is an example of initial curves for a synthetic image, andFIG. 64 b is an example of corrected curves for a synthetic image;

FIGS. 65 a, 65 b and 65 c show an example of shape recovery of a curvewhen the external forces are removed, and FIG. 65 d is an example ofboth initial curve (in white) and final curve (in black);

FIG. 66 is a schematic flow chart showing how to build an illustrativeembodiment of the computational image model according to the invention;and

FIG. 67 is a schematic flow chart showing how to build a computationalframework for solving a problem, in accordance with an illustrativeembodiment of the present invention and using the computational imagemodel of FIG. 66.

DETAILED DESCRIPTION OF THE ILLUSTRATIVE EMBODIMENTS

The illustrative embodiments of the present invention are generallybased on a data structure for images based on n-pixels, in which theimage support, the image quantities and the allowable operations arespecified separately. In this data structure, mathematics and physicsare unified; that is, the data structure allows taking into accountconstraints originating in physics, mathematics, and the future use ofthe image. The image dimension is explicit, which allows us to designalgorithms that operate in any dimension. The foregoing and otheradvantages and new open problems of this data structure will bediscussed herein below.

One of the goals of the present invention is to provide a computationalimage model or a data structure that is capable of capturing all objectproperties that are needed for a given task. FIG. 66 is a schematic flowchart summarizing the procedure (successive steps 6601-6606) forbuilding such a computational image model according to an illustrativeembodiment of the invention. A data structure of the image is the formalspecification of the image variables, image quantities, the associationbetween image quantities and variables that enables capture of thegeometrical and topological properties of objects as well as theirphysical and mathematical behavior. The abstract view of an image as adata structure as is used in computer programs is defined by theattributes and a collection of meaningful operations. The attributes arethe image support and quantities that are assigned to the image supportsuch as the image radiometry (e.g., color and grey level) or any featurethat can be deduced from the radiometry (e.g. texture). These quantitiesare scalar, vector or tensor. The allowable operations are of two kinds:the operations that are problem independent such as read and write andthose that are problem dependant such as object deformation, diffusionand optical flow. To summarize, the image is defined by the support,quantities and allowable operations. The latter three iterns will bedefined in detail in the following description.

Image Support

Often, discretization of an image is achieved via a cubic tessellation.For example, a 2D (two-dimensional) image support is formed by unitsquares commonly called pixels. Similarly, a 3D (three-dimensional)image is a tessellation of unit cubes commonly called voxels. Moregenerally, the illustrative embodiments of the present invention willconsider the image in n dimensions as a set of unit n-cubes, which arecommonly called n-pixels. When n=0, the image is a set of points; whenn=1, a set of edges; when n=2, a set of squares; when n=3 a set ofcubes, and so on. Any two n-pixels are either disjoint or intersectthrough a common i-pixel, where i<n. This subdivision of the imagesupport is not unique. Several other geometrical forms such as, forexample, triangles or hexagons can be used. It has been proven that thetopological features of the image support do not depend on thesubdivision used [13]. The cubical subdivision is commonly used in imageprocessing and computer vision and will therefore be used in thenon-restrictive illustrative embodiments of the present invention. Onedoes not need to explicitly define an orientation of pixels. Indeed, thedefinition of a pixel as a product of intervals in a certain ordercoupled with the natural order of real numbers imposes an orientation oneach coordinate axis and also on the canonical basis of R^(n) (see FIG.1).

Formally, a pixel σ ε R^(n) as a geometric entity is translated into analgebraic structure as follows:σ=I ₁ ×I ₂ × . . . ×I _(n)   (1)where x is the Cartesian product and I_(i) is either a singleton or aninterval of unit length with integer endpoints; i.e., I_(i) is eitherthe singleton {k}, in which case it is said to be a degenerate interval,or the closed interval [k,k+1] for some integer k. The number q ε {0,1,. . . ,n} of non-degenerate intervals in Equation (1) is, by definition,the dimension of σ, which will be referred to as a q-pixel. For q≧1, letj={k₀,k₁, . . . ,k_(q−1)} be the ordered subset {1,2, . . . ,n} ofindices for which I_(k) _(j) =[a_(j),b_(j)] is a non-degenerateinterval. DefineA _(k) _(j) σ=I ₁ × . . . ×I _(k) _(j) ⁻¹ ×{a _(j) }×I _(k) _(j) ₊₁ × .. . ×I _(n)andB _(k) _(j) σ=I ₁ × . . . ×I _(k) _(j) ⁻¹ ×{b _(j) }×I _(k) _(j) ₊₁ × .. . ×I _(n)where A_(k) _(j) σ and B_(k) _(j) σ are, respectively, the front(q−1)-face and the back (q−1)-face of σ. Each of these (q−1)-faces is a(q−1)-pixels. These faces are then called (q−1)-faces of σ. In the samemanner, one can define the (q−2)-faces, . . . , down to the 0-faces ofσ. FIG. 1 shows a 2-pixel A, with its 1-faces a, b, c, d. The 0-faces ofA are the vertices that are not represented here for the sake of clarityof the picture. The boundary of the q-pixels σ enables the writing ofthe relationship between a q-pixels and its (q−1)-faces in algebraicform. By definition this is the alternating sum of its (q−1)-faces, i.e.$\begin{matrix}{{\partial_{q}\sigma} = {\sum\limits_{i = 0}^{q - 1}{\left( {- 1} \right)^{i}\quad\left( {{B_{k_{i}}\sigma} - {A_{k_{i}}\sigma}} \right)}}} & (2)\end{matrix}$

For example, the boundary of A in FIG. 1 is given by the followingrelation:∂₂σ=(−0×[0,1]+1×[0,1])−([0,1]×0−[0,1]×1)=(c−a)+(d−b).

So far, the pixel, its faces, and the association between them have beendefined geometrically and algebraically. Now, the image support will bedefined as a geometrical entity, called a cubical complex, and then itsalgebraic structure will be written. A cubical complex K in R^(n) is acollection of q-pixels where 0≦q≦n such that:

-   -   Every face of a pixel in K is also located in K; and    -   The intersection of any pair of two pixels of K is either empty        or formed by a common face of both pixels of the pair.

The first condition implies that all faces of a pixel belong to thecubical complex.

The second condition concerns the organization of the cubical complex.If the intersection of two pixels is empty, then the image support isformed by one or more connected components. This condition provides animage support that is more general than existing definitions since itallows a formal specification of a cubical complex that is formed eitherby one or more image supports (e.g., an image sequence) or by severaldistinct binary objects. When the intersection is a face, certaingeometrical configurations of the complex are ruled out. For example,FIG. 2 a illustrate a subdivision that is a cubical complex and FIG. 2 billustrates a subdivision that is not a cubical complex.

The dimension of the cubical complex K is, by definition, the largestnumber q for which K contains a q-pixel.

As in the case of the q-pixel, the cubical complex can be written inalgebraic form. Given a topological space X ⊂ R^(n) in terms of acubical complex, the set of all q-pixels of X is denoted E_(q)={σ_(q) ¹,. . . ,σ_(q) ^(N) ^(q) }. A q-chain in X is a formal sum of integermultiples of elements of E_(q). More precisely, it is a linearcombination $\begin{matrix}{{\lambda_{1}\sigma_{q}^{1}} + \cdots + {\lambda_{N_{q}}\sigma_{q}^{N_{q}}}} & (3)\end{matrix}$where λ₁, . . . ,λ_(N) _(q) are integers. For example, in FIG. 1,a−c+d−b is a 1-chain. Two q-chains are added by adding correspondingcoefficients. The set of q-chains can be given the structure of a freeAbelian group with basis E_(q), usually denoted by C_(q)(X). Since onlyfinite complexes will be considered, the groups C_(q)(X) are finitelygenerated and C_(q)(X)=0 if q is greater than the dimension of thecomplex; naturally, C_(q)(X)=0 if q<0.

To define the chains that are associated With the faces of pixels of acubical complex, Equation (2) is extended by linearity to all q-chainsto obtain the boundary map ∂_(q):C_(q)(X)→C_(q−1)(X). Note that ∂₀=0since C⁻¹(X)=0. The boundary map satisfies the following property [9]:∂_(q) ∘ ∂_(q+1)=0   (4)

To summarize, the discrete image support of any dimension n is formed byn-pixels. Unlike conventional image models, the n-pixel is a dimensionalgeometric entity formed by other geometric entities called faces. Thegeometrical n-pixel is translated into a recurrent algebraic structure,more reliable for mathematical handling. All n-pixels of the imagesupport form a cubical complex, a geometrical entity that is translatedinto an algebraic structure called the chain. The association betweenthe n-pixel and its faces, and hence between chains of successivedimensions is given by a boundary operator. It should be noted that theuse of any other geometrical primitive such as a triangle or hexagon forsubdividing the image support affects the computational complexity ofthe derived algorithms, but it affects neither the topological featuresof the image support nor the computation rules for these features.

Image Quantities

In the foregoing description, the concept of the finite cubical complexwas introduced to give an algebraic description of the discrete imagesupport. A similar formulation is needed to describe the image field(scalar, vector, matrix) over the discrete image support. For thispurpose, let's return to the above described chain concept to give amore general definition. Considering a cubical complex K of dimension n,each q-pixel (q≦n) of K is associated to a coefficient in the ring(A,+,*), where the elements may be scalars (gray level), vectors (color,gradient), matrices (Hessian), etc. The chain is the formal sum$\sum\limits_{i}{\lambda_{i}c_{q}^{i}}$where λ_(i) is a coefficient in (A,+,*), and the generators c_(q) ^(i),∀i form a basis of an Abelian group. The chain can be seen as a vector[λ₁, . . . ,λ_(N) _(q) ]^(T), where N_(q) is the number of generatorsused. Consequently, two chains can be summed by adding theircoefficients and multiplied using the scalar product. The addition andmultiplication are taken according to the definition of a ring.Moreover, one cam define a null chain (respectively unit chain) whosecoefficients are all equal to the null (respectively unit) element ofthe + (respectively *) operation of (A,+,*). Consequently, q-chainsdefine an image model that has attractive computational properties sincethey form a rich algebraic structure, the module (i.e, a vector spacedefined on a ring).

To briefly show how to use the chain model in image processing, a simpleillustrative example concerning any global transform of the image suchas histogram equalization or thresholding is presented. The chaincoefficients are the gray levels. The formal expression of the globaltransform implies two applications: {overscore (H)}:(ε_(A),+,*)→(ε_(A),+,*) and H: (A,+,*)→(A,+,*), where (ε_(A),+,*) is amodule. They are defined by${\overset{\_}{H}\left( {\sum\limits_{i}{\lambda_{i}c_{i}}} \right)} = {\sum\limits_{i}{{H\left( \lambda_{i} \right)}{c_{i}.}}}$

The drawback of chain-based image models is that the physical ormathematical quantities and the image support are described together ina formal sum. Consequently, the chain coefficients combine mathematicalor physical quantities, pixel orientation, and possibly other quantitiessuch as weights associated with pixels. For example, there is ambiguityin the interpretation of the sign of λ_(i): it may be due to theorientation, the physical quantities, the weights or theirmultiplication. This image model can be considered adequate, especiallyin physics [14], engineering and computer graphics [11, 5], where chainshave been used to model quantities. However, this illustrativeembodiment of the present invention proposes to refine this quantitymodel to overcome the confusion produced by the chain coefficients.

To reflect only the geometry of the image support (e.g., orientation andmultiplicity), what follows will consider q-chains with integercoefficients. We are looking for an application F_(q): C_(q)(X)→(B,+,*),which associates a global quantity (energy, gray level, color, flux,tensor, etc.) with all q-pixels, where q≦n and (B,+,*) is a ring. As inthe case of the chain-based image model, for two adjacent q-pixels c_(q)¹ and c_(q) ², F_(q) must satisfy F_(q)(λ₁c_(q) ¹+λ₂c_(q)²)=λ₁F_(q)(c_(q) ¹)+λ₂F_(q)(c_(q) ²), which means that the sum of thequantities generated within each q-pixel is equal to the quantitiesgenerated within the two q-pixels. F_(q) can be extended by linearity toany q-chain ${\sum\limits_{i}{\lambda_{i}c_{q}^{i}}},$where each λ_(i) is an integer, as follows:${F_{q}\left( {\sum\limits_{i}{\lambda_{i}c_{q}^{i}}} \right)} = {\sum\limits_{i}{\lambda\quad{{F_{q}\left( c_{q}^{i} \right)}.}}}$

F_(q) is called a q-cochain. To illustrate the cochain concept, let usconsider a 2-pixel c₂ and a vector field V. A 0-cochain is defined bythe value of V at 0-pixels. A 1-cochain is ∫_(∂c₂)  V ⋅ 𝕕s,the line integral over the faces of c₂. A 2-cochain is ∫_(c₂)  V ⋅ 𝕕S,the surface integral over the 2-pixel c₂, and the “.” the dot product.

To summarize, a cochain allows us to associate quantities with theq-pixels and with the faces thereof. Unlike existing image models, themodel according to the illustrative embodiments of the present inventionprovides a rich structure that allows the definition of both local andglobal quantities.

Operations

A generic operation that can be instantiated depending on the problemthat we are dealing with will now be defined. Having in mind that aq-pixel has 3^(q) faces, the generic operation should specify thealgebraic relationship between the quantities (i.e., cochains)associated with these faces. Based on the linearity principle, thequantity of a given q-pixel is transferred to its cofaces with the sameor opposite sign, according to the agreement of its orientation and theorientation of its cofaces. The quantities that are transferred to the(q+1)-pixels by its faces are summed. More formally, it should bereminded that the relationship between the (q−1)-chain and the q-chainis given by the boundary operator. Similarly, the relationship betweenthe q-cochain and the (q+1)-cochain is given by the coboundary operatorδ_(q): C^(q)→C^(q+1), where C^(q) is the Abelian group of q-cochains.Given a (q+1)-chain c, this coboundary operator is defined by:δ_(q) F _(q)(c)=F _(q)(∂_(q+1) c)   (5)

The capacity of the cochain and the coboundary to model a given problemwill now be discussed. The cochain is a linear application and shouldfulfill the coboundary requirement of Equation (5). Thus a questionconcerns the limits in modeling a given quantity. It is difficult toanswer this question for the general case. Much investigation is neededfirst. In the present illustrative embodiment, this question is onlyanswered by identifying several problems that can be modeled by thecochain and coboundary. Certain problems such as convolution and itsapplications (smoothing, numerical differentiation, high-pass filtering,noise estimation) and the Fourier transform can be modeled by thecochain without approximation since they only require setting thecoefficients λ_(i) of the cochain to specific values (see [16] for thecase of numerical differentiation). Thresholding can be represented bythe cochain without approximation since${{H\left( {\sum\limits_{i}{\lambda_{i}Q_{i}}} \right)} = {\sum\limits_{i}{\lambda_{i}{H\left( Q_{i} \right)}}}},$where H: (B,+,*)→(B,+,*). Other problems can be broken down to basiclaws, each of which can be described by the topological Equation (5).For example, it has been pointed out [14] that many physics problems canbe broken down into basic physical laws such as balance and constitutivelaws. As it will be showed herein below, balance law can be written in adiscrete form by using the topological Equation (5) withoutapproximation. The constitutive laws cannot be translated in algebraicform without approximation. Usually, they require the link betweencochains that belong to different cubical complexes. For example, in thecase of dual complexes, two cochains are linked by an algebraic linearsystem. This transformation is called “codual operation”. Moregenerally, 0-cochains represent local quantities. However, q-cochains(q>1) represent global quantities (e.g., an integral or the summation ofa differential form) since they are associated with the algebraicstructure of an edge, an area, a volume, etc. Thus, the cochain,coboundary, and codual are generic algebraic structures that can beinstantiated by physical or mathematical laws. The exact translation ofgiven problem in terms of q-cochains is possible if one is able to findthe basic laws that can be described without approximation by eithercochains or the topological Equation (5).

In the previous example, the coboundary can be interpreted as follows:assuming that the vector field is conservative $\begin{matrix}{{V = {\Delta\quad v}},} & {{\int_{c_{1}}^{\quad}{\Delta\quad{v \cdot {\mathbb{d}s}}}} = {{v(b)} - {v(a)}}}\end{matrix}$constitutes a coboundary operator since it may be written asδ₀F₀(c₁)=F₀(∂₁c₁), where a and b are the faces of c₁. Similarly,∫_(c₂)  div(V)𝕕S = ∫_(∂c₂)  V ⋅ n  𝕕s,where n is the outward unit normal vector to ∂c₂, constitutes acoboundary operator since it may be written as δ₁F₁(c₂)=F₁(∂₂c₂). Thisexample shows that the coboundary operator may be an exact discreterepresentation of the fundamental calculus theorems (line integral andGauss).

Thanks to the chain, cochain, coboundary and codual concepts, the imagemodel of the illustrative embodiments of the present invention can takeinto account the mathematical or physical laws related to theapplication. It can thus be instantiated to the various problems ofimage processing and computer vision. The underlying computationalframework is strongly different form existing frameworks. For example,let us consider physics based problems such as optical flow, diffusionand deformation. Existing frameworks can be summarized as follow: 1)modeling by partial differential equation; 2) resolving the PDE by usingnumerical analysis scheme or Fourier space.

The computational framework, according to the illustrative embodimentsof the present invention, can be summarized as follow: 1) identificationof basic laws associated to the problem (Block 6701 of FIG. 67); 2)definition of an image support including the number of cubical complexesand the dimension of the cubes (e.g., for the case of a multi-resolutionprocessing) (Block 6702 of FIG. 67); 3) definition of global and localquantities (Block 6703 of FIG. 67); 5) instantiation of the coboundaryand codual operators (Block 6704 of FIG. 67); and 6) the resolution ofresulting algebraic system (Block 6705 of FIG. 67). The advantages ofthis framework are described hereinbelow and will be better defined intwo practical examples.

To summarize, the coboundary operator links quantities associated withthe faces of an n-pixel. Codual operator links quantities associated tocomplexes of an image. If a given problem can be broken down into basiclaws, the cochain and coboundary are the discrete representation ofthese basic laws. Cochain, coboundary, and codual are generic conceptsthat can be instantiated by physical or mathematical laws. Thus, theycan be used in various computer vision and image processing problems.

AN ILLUSTRATIVE EXAMPLE

Let us consider the linear isotropic diffusion in gray level imageswhich is a physics-based problem. One can find all details in reference[4]. For the sake of clarity and without loss of generality, theanalysis will be limited to considering the 2D global differentialequation for heat flow in a homogeneous medium. Let V be a 3D regionwith boundary S inside the flow. The rate of heat flow across boundary Sout of V is given by: $\begin{matrix}{{\int{\int{\int_{V}^{\quad}{{\sigma\left( {x,t} \right)}\quad{\mathbb{d}V}}}}} = {- {\int{\int{\int_{V}^{\quad}{{\nabla{\cdot \left( {\lambda\quad{\nabla{\mathcal{T}\left( {x,t} \right)}}} \right)}}\quad{\mathbb{d}V}}}}}}} & (6)\end{matrix}$where x is a spatial vector, t represents time, and λ is a positiveconstant.

To resolve this problem, the image support is defined by a continuousscalar field T, the temperature (i.e., gray level), two dual cubicalcomplexes (i.e., two chains), and three cochains. If only one cubicalcomplex is used, two different orientations are associated with each1-face. To overcome this problem, two dual complexes (primary andsecondary) are used (see FIG. 3). Concerning the use of three cochains,it has been pointed out in reference [4] that this EDP is formed bythree basic physical laws. Each cochain is associated with one law. Thefirst is the thermal tension law (also known as Fick's Law), whichstates that heat flows from regions of higher temperature towardsregions of lower temperature. The direction of the gradient ∇T is thedirection of the largest increase in temperature; the heat flows in theopposite direction. Formally, this law is written as follows:g(x, t)=−∇T(x, t)   (7)

The primary cubical complex is the support for this balance law. Theorientation plotted on this cubical complex is the direction of the pathon which the integral is computed. Let us assume that the 0-cochain isthe temperature T associated with 0-pixel c₀. A 1-cochain G associatedwith 1-pixel c₁ is the global thermal tension transferred by the two0-pixels that are the faces of c₁. Consequently, the topologicalequation is:G(c ₁)=δ₀ T(c ₁)=T(∂₁ c ₁)   (8)

By the linearity of cochains, this topological equation is valid for all1-chains.

The second law, called the heat source law, concerns the net outflow ofinternal thermal energy at the point x and time t. This is a balancelaw. It is given by:σ(x, t)=∇.q(x, t)   (9)

When ∇.q(x,t)>0, the outflow is positive and thermal energy must flowaway from x. Similarly, if ∇.q(x,t)<0, the inflow is larger than theoutflow and thermal energy increases at x. The equilibrium for adiffusion process is attained if ∇.q(x,t)=0.

Let us consider the secondary cubical complex. The orientation plottedon this cubical complex is the direction of the flow. The 2-cochain Σassociated with the 2-pixel c₂ is the global heat transferred by thefaces of c₂Σ(c ₂)=δ₁ Q(c ₂)=Q(∂₂ c ₂)   (10)

This topological equation is valid for all 2-chains. The third law isconstitutive (it depends on the medium feature). It makes the linkbetween the flow density law and the heat source law. It is given by:q(x, t)=λg(x, t)   (11)

This equation cannot be discretized exactly, since it involves globalquantities defined on two dual complexes. Consequently, the 2-cochain±cannot be computed without approximation from the 1-cochain G. Forexample, they can be linked as a linear equation system ΛQ=±, where Λ isthe coefficient matrix. FIG. 4 a gives an example of original image andFIG. 4 b and example of image smoothed using this computational scheme.

The data structure associated to the linear diffusion problem is definedby: 1) two dual cubical complexes; 2) two cochains G and Λ for globalquantities and a scalar field T; 3) two coboundary operations forbalance laws and a codual operation that represents the constitutivelaw. The framework is summarized as follows: g is approximated by abilinear polynomial. The cochain G is computed using a line integral. qis computed using the constitutive law in Equation (11). The cochain ±is computed using Gauss's theorem. Finally, a system of linear equationsobtained from Equation (11) is obtained where the unknown variables areT. It should be noted that the cochains in Equations (8) and (10) arecomputed without approximation.

The image model according to the illustrative embodiments of the presentinvention and described hereinabove may generally be characterized bythree major points: 1) the image support and quantities are definedseparately and then linked together via algebraic language; 2) the pixelis dimensional and is written in an algebraic form; 3) both local andglobal quantities are represented by the cochains and coboundaryoperators.

Each of these specificities will now be discussed to show theirstraightforward consequences for image processing.

The separability of the image model allows a distinction between imagevariables and image quantities. The image variables offer numerouspossibilities for existing mathematical formulations such as the use ofalgebraic topology to help in the design of algorithms. For example,binary image algorithms are written as algebraic systems [16]. Thewell-defined quantities allow the use of physics, vector analysis ordifferential forms in the design of algorithms. Taking image support andimage quantities together, well-known problems such as those ofdiffusion and optical flow in gray level images can be written asalgebraic systems. Furthermore, the transfer of quantities between agiven domain and its boundary is straightforward, using the concepts ofcochain and coboundary as a general framework. For example, as shownhereinabove, in vector calculus, this transfer is easier thanks to thethree fundamental calculus theorems, the line integral, Stokes'stheorem, and Gauss's theorem.

Unlike existing image models, by considering the pixel as dimensionalprimitive the connectivity paradox of the image support is avoided [8].That, is, the well-known Jordan theorem is fulfilled. Its decompositioninto faces and the use of cubical structures such as chains make thedimension of the image explicit. Algorithms designed according to thisformalism operate in any dimension. This fact overcomes the traditionallimitations that are faced in designing an algorithm, say in onedimension, extending it to two dimensions, and then to three dimensions,and so on. Each extension step may be a difficult task.

The definition of the cochain depends on the problem that is dealt with.Thus, this image model offers real flexibility for the integration ofmathematical objects (scalar, vector, tensor) and physical laws(balance, constitutive). Furthermore, the use of global quantitiesassociated with an n-pixel implies noise reduction. In fact, globalquantities are computed from the field by using the integral or thediscrete summation. As the opposite operation from the differentiation,which enhances high frequencies of the image, the integral performs asmoothing operation. It allows us to reduce the order of the derivativeused in an image-processing scheme. Consequently, the introduction ofglobal quantities may allow the use of higher-order derivatives.

Another contribution of the image model according to the illustrativeembodiments of the present invention concerns the numerical scheme usedto solve nonlinear problems such as the diffusion problem and elasticmatching. It should be reminded that, usually, a problem is formulatedand then a numerical analysis scheme is used. The numerical analysisscheme may not have been derived for exactly this formulation and manyapproximations must be made. The explanation of intermediate results isnot available. Consequently, no clear idea is available about theconvergence of the numerical analysis scheme and the numerical resultsobtained may be a broad approximation of the desired solution. Based onthe problems tackled, it is concluded that in the image model presentedhere, the numerical scheme is deduced from the problem model with-littleor no approximation [4, 12]. In fact, various problems may be brokendown into basic laws and then reformulated in terms of cochains andcoboundaries. Thus they can be written as linear algebraic systems andsolved.

PRACTICAL EXAMPLE #1 Physics-Based Resolution of Diffusion and OpticalFlow

An alternative to Partial Differential Equations (PDES) will now bedescribed for the solution of three problems in computer vision: linearisotropic diffusion, optical flow and nonlinear diffusion. These threeproblems are modeled using the heat transfer problem. Traditionally, themethod for solving physics-based problems such as heat transfer is todiscretize and solve a PDE by a purely mathematical process. Instead ofusing the PDE, the global heat problem can be decomposed into basiclaws. It will be demonstrated that some of the basic laws admit an exactglobal version since they arise from conservation principles. It willalso be showed that the assumptions made on the other basic laws can bemade wisely, taking into consideration knowledge about the problem andthe domain. The above-described image model will be used to allowencoding of physical laws by linking a global value on a domain withvalues on its boundary. The resulting algorithm performs in anydimension. The numerical scheme is derived in a straightforward way fromthe problem modelled. It thus provides a physical explanation of eachsolving step in the solution.

Background

In recent years, Partial Differential Equations (PDEs) have attractedincreasing interest in the field of computer vision. Since PDEs havebeen the subject of much study by numericians, powerful numericalschemes have been developed to solve them. Consequently, domains such asimage enhancement, restoration, multi-scale analysis and surfaceevolution all benefit greatly from PDEs [25].

One important class of equations governing certain physical processes isthe linear elliptic PDE of the general form known as the Helmholzequation:∇² u(x)+p(x)u(x)=f(x)   (12)where x denotes a vector in the n-dimensional space, u(x) is thedependent variable, ∇² is the Laplacian operator, and p(x) and f(x) arespatial functions. When p(x)=0, this corresponds to the Poisson equation[21, 32] (also known as the non-homogeneous Laplace equation [30, 44]).One of the physical processes governed by Equation (12) is steady-stateheat transfer.

In the field of computer vision, Equation (12) may arise from twoapproaches. The first is variational calculus. As a matter of fact, manyproblems such as shape from shading [39], surface reconstruction [32,40] and the computation of optical flow [29] can be formulated asvariational problems. The solutions to these variational problems aregiven by Euler-Lagrange equations, which are in the form of Equation(12) [31, 32]. The second approach is physics-based. For example,diffusion processes arise from heat equations and shock filters fromwork in fluid mechanics [25]. For both the variational and thephysics-based approaches, the resulting PDEs are continuous and have tobe discretized.

Traditionally, the discretization of PDEs in computer vision has beendone by applying finite difference methods [23, 31, 39, 40]. Equation(12) is solved iteratively using either a direct Fourier-based Poissonsolver for each iteration [39], finite elements [24], or spectralmethods [32]. Iterative methods such as those in [39] do not ensureconvergence unless smoothness is very high [21].

Existing methods for the resolution of problems involving PDEs can besummarized as follows: 1) identification of basic laws; 2) combinationof the basic laws in order to write the PDEs; 3) discretization of thePDEs; 4) resolution of the PDEs via a numerical method.

This process, which has been used in various fields of application, ispurely mathematical. Consequently, it has the following drawbacks: 1)Some quantities involved in the solution process do not have a physicalinterpretation; 2) This lack of interpretation is manifested inintermediate solutions involving iterative processes and since thesesolutions cannot be physically explained, discovery of the optimalsolution cannot be ensured in an optimal time.

Solution

To overcome these drawbacks, an alternative to PDE resolution in thecontext of the heat transfer problem is proposed and will be describedhereinbelow.

Generally; basic laws in physics-based problems are combined into aglobal conservation equation [42] that is valid on-the whole body or apart of it. A local conservation equation (PDE) is then obtained byconsidering the particular case of a part of a body reduced to a point.

In discrete problems such as those encountered in computer vision, thecontinuous domain is subdivided into many sub-domains in which there isonly one value available, which can be considered as a global value.Therefore, instead of using the PDE, this illustrative embodiment of thepresent invention proposes to use the global conservation equationdirectly on each sub-domain.

In order to handle these physical laws which link global values atpoints, lines, surfaces, volumes, etc. the image model with roots incomputational algebraic topology described hereinabove is used. Thismodel makes it possible to represent global values on entities of anydimension at the same time.

The above described methodology presents a number of advantages:

-   -   1) Many of the basic laws arise out of conservation principles        and hence they are valid either at a point (local form as in        Equation (12)) or over an entire region (global form).        Fundamental theorems of calculus such as the Gauss, Green and        Line Integral theorems allow the computation of the coboundary        operator without any approximation.    -   2) Some laws require approximations that can be performed        wisely, taking into account knowledge about the problem and the        domain.    -   3) The intermediate results have a physical explanation because        they represent physical quantities. For that reason, every step        has a physical interpretation. Thus there are no longer problems        of non-optimality of the solution, because we avoid non-temporal        iterative methods.    -   4) As mentioned earlier, this method can be used with other        physics based problems by applying the appropriate basic laws.    -   5) Thanks to the image model, the resulting algorithm performs        in any dimension.    -   6) The computational rules associated with the coboundary        operator can be changed without changing the formalism of the        operation itself.    -   7) The same formalism can be used for pixel-based and other        types of decomposition of the image (e.g. regions).

In order to validate the method, the equation for steady state heattransfer is resolved in two applications: the linear diffusion andoptical flow. These problems generate equations of the form of Equation(12) or its global version. The present methodology can also be used toresolve unsteady heat transfer with no source and we apply it tonon-linear diffusion.

Physical Principles (Explanation of the Physical Foundations of the HeatTransfer Problem)

Two interesting particular cases for diffusion and optical flow problemscan be distinguished: steady-state heat transfer and unsteady heattransfer with no source.

Two important classes of laws are present: conservation and constitutivelaws. Conservation laws are independent of the properties of thematerial, whereas constitutive laws are specific to them. The physicalproperties associated with a moving object are energy, work and heat. Inwhat follows, each of these quantities are described.

Energy Modeling

Some quantities for a continuous body occupying a volume V bordered by asurface S in a 3D space will first be defined. Such a body can be saidas composed of an infinite number of particles (as many as desired),these particles being the smallest elements [33]. FIG. 5 a illustratessuch a body. At time t, a particle labelled X occupies a specificposition:x=(x(t), y(t), z(t))

Each particle can move in space, so a velocity vector is associated withit at time t:${v^{*}\left( {X,t} \right)} = {\frac{\mathbb{d}x}{\mathbb{d}t} = {v\left( {x,t} \right)}}$

Physical quantities can be associated with a particle labelled X(material description) or a position x in space (spatial description).For example, v*(X,t) is the material description of the velocity ofparticle X and v(x,t) is the spatial description of the particle locatedat position x. For the present purpose, spatial descriptions are used toderive the heat transfer equation. Vector n(x,t) is the outwarddirection of the surface at point x.

The mass Δm of a small amount of volume ΔV of a body is a measure of itsinertia (tendency to resist motion). The term mass density, ρ is used todenote the following quantity:$\rho = {\lim\limits_{{\Delta\quad V}\rightarrow 0}\frac{\Delta\quad m}{\Delta\quad V}}$ρ(x,t) is thus the mass density of the particle located at x at time t.

Two kinds of energy are associated with a moving object: kinetic andinternal energy. Kinetic energy is a measure of the state of motion of abody: the faster the body moves, the greater its kinetic energy [28].Because it is a measure of inertia, kinetic energy also fakes the massinto account. For a particle located at x at time t, the kinetic energyis thus defined as${K\left( {x,t} \right)} = {\frac{1}{2}{\rho\left( {x,t} \right)}\quad\left( {{v\left( {x,t} \right)} \cdot {v\left( {x,t} \right)}} \right)}$where “·” is the dot product. To obtain the kinetic energy for theentire body at time t, K(x,t) is integrated over the volume V:${K(t)} = {\frac{1}{2}{\int{\int{\int_{V}^{\quad}{{\rho\left( {x,t} \right)}\quad\left( {{v\left( {x,t} \right)} \cdot {v\left( {x,t} \right)}} \right)\quad{\mathbb{d}V}}}}}}$where dV is an infinitesimal amount of the volume V.

Internal energy is a measure of the state of temperature of a body. Thehotter the body, the greater its internal energy. At time t, eachparticle has an internal energy density ε(x,t) associated with it. Theinternal energy density is proportional to the temperature of theparticle T(x,t) with a material-specific heat constant c; that is,ε(x,t)=cT(x,t). For the entire body, the total internal energy isintegrated over the volume V: $\begin{matrix}{{E(t)} = {\int{\int{\int_{V}^{\quad}{{\rho\left( {x,t} \right)}\quad{ɛ\left( {x,t} \right)}\quad{\mathbb{d}V}}}}}} & (13)\end{matrix}$

The total energy for the body can now be defined as K(t)+E(t).

Work Modeling

Let us suppose that a body is submitted to an external force f_(e)(x,t)(e.g. a traction force) and an internal density force b(x,t) (e.g.gravity). FIG. 5 b presents the action of external and internal forces.Work is defined as the energy transferred to a body by means of a forceacting on the body [28]. Work is negative when the energy is transferredfrom the body. Suppose that F(x) is an internal or external force thatis constant over time, acting on a particle located at x during anamount of time t. This force will produce a displacement of the particleto position x₁. This displacement is Δx=x₁−x. The work w(F, x) done bythis force during this time is:w(F, x)=F(x)·Δxand the instantaneous power P(x, t) is: $\begin{matrix}{{P\left( {x,t} \right)} = {\lim\limits_{{\Delta\quad t}\rightarrow 0}\frac{w\left( {F,x} \right)}{\Delta\quad t}}} \\{= {{F(x)} \cdot {\lim\limits_{{\Delta\quad t}\rightarrow 0}\frac{\Delta\quad x}{\Delta\quad t}}}} \\{= {{F(x)} \cdot \frac{\mathbb{d}x}{\mathbb{d}t}}} \\{= {{F(x)} \cdot {v\left( {x,t} \right)}}}\end{matrix}$

External forces act essentially on the surface of the body. Theinstantaneous work P_(e)(t) done by the external forces on the entirebody is thus the result of the integration of the external power overthe surface: P_(e)(t) = ∫∫_(S)  f_(e)(x, t) ⋅ v(x, t)  𝕕Swhere dS is an infinitesimal part of the surface of the body.

Defining b(x, t) as the internal density force, the internal force isthus ρ(x, t)b(x, t). The rate of work over time done by the internalforces on the entire body is obtained by integrating internal work overthe volume: P_(i)(t) = ∫∫∫_(V)  ρ(x, t)  (b(x, t) ⋅ v(x, t))  𝕕VThe total work is thus P(t)=P_(e)(t)+P_(j)(t)Heat Modeling

Heat can be defined as the energy transferred to a body owing to adifference in temperature. The heat flow density vector q(x, t) is ameasure of the rate of heat conducted into the body per unit area perunit of time. How q(x, t) is defined will be explained later. Theexternal heat addition rate over time is the amount of heat coming fromoutside the body and entering by its surface. It is computed byprojecting q(x, t) onto the inward normal vector (−n(x, t)) andintegrating this projection over the surface: $\begin{matrix}{{Q_{e}(t)} = {\int{\int_{S}^{\quad}{{{q\left( {x,t} \right)} \cdot \left( {- {n\left( {x,t} \right)}} \right)}\quad{\mathbb{d}S}}}}} & (14)\end{matrix}$

Now if a body has a rate of heat generation per unit of volume and timer(x, t), the internal rate of heat addition over time is computed byintegrating r(x, t) over the volume:Q_(i)(t) = ∫∫∫_(V)  ρ(x, t)  r(x, t)  𝕕V

FIG. 5 c shows q(x, t) and r(x, t). To simplify, the source σ(x, t) isdefined as the rate of heat generated in a particle located at x perunit of volume and time:σ(x, t)=ρ(x, t)r(x, t)   (15)

In many cases, this source is known. However, it could also be a linearfunction of the temperature: σ(x, t)=a(x, t)+b(x, t)T(x, t) [38]. Thetotal rate of heat addition over time is thus:Q(t)=Q _(e)(t)+Q _(i)(t)Energy Conservation Law

A class of equations in continuum mechanics are those describing theconservation (equilibrium) principles. They express the conservation ofcertain physical quantities (mass, momentum, energy, etc.) over anentire body and as such they take the form of global equations over thewhole body or a part of it [33].

Conservation principles can be seen intuitively as follows: the changein the total amount of a physical quantity inside a body is equal to theamount of this quantity entering or leaving the body (through theboundary) and the amount generated or absorbed within the body. Theselaws are applicable for all continuous materials, moving and stationary,deformable and non-deformable, and must always be satisfied. The globalconservation equations can then be used to derive their localcounterparts, called the associated field equations, which are valid ateach point of the body including its borders.

The first law of thermodynamics, which is relevant for the understandingof the heat transfer equation will now be discussed. This law involvesboth kinetic and internal energies and states that the total variationof energy in a body (or a part of a body) is the result of the time rateof work and the rate of heat addition combined: $\begin{matrix}{{\frac{\mathbb{d}\quad}{\mathbb{d}t}\left( {{E(t)} + {K(t)}} \right)} = {{P(t)} + {{Q(t)}.}}} & (16)\end{matrix}$

For heat transfer, the only interest resides in the case of immobilebodies, that is v=(0,0,0), n(x,t)=n(x) and ρ(x,t)=ρ(x). Equation (16)thus becomes:${\frac{\mathbb{d}\quad}{\mathbb{d}t}{\int{\int{\int_{V}^{\quad}{{\rho(x)}\quad{ɛ\left( {x,t} \right)}\quad{\mathbb{d}V}}}}}} = {{\int{\int_{S}^{\quad}{{- \left( {{q\left( {x,t} \right)} \cdot {n\left( {x,t} \right)}} \right)}\quad{\mathbb{d}S}}}} + {\int{\int{\int_{V}^{\quad}{{\sigma\left( {x,t} \right)}\quad{\mathbb{d}V}}}}}}$which now states that the thermal energy variation in a body is due tointernal heat production added to the heat flowing into the body. Usingthe divergence theorem for Q_(e) [33] and recalling that ε(x,t)=cT(x,t),we obtain the thermal energy conservation law: $\begin{matrix}{{\int{\int{\int_{V}^{\quad}{{\rho(x)}\quad c\quad\frac{\mathbb{d}\quad}{\mathbb{d}t}{T\left( {x,t} \right)}\quad{\mathbb{d}V}}}}} = {{\int{\int{\int_{V}^{\quad}{{{- \nabla} \cdot {q\left( {x,t} \right)}}\quad{\mathbb{d}V}}}}} + {\int{\int{\int_{V}^{\quad}{{\sigma\left( {x,t} \right)}\quad{\mathbb{d}V}}}}}}} & (17)\end{matrix}$where ∇. is the divergence operator. To simplify, let us define thetemperature variation${h\left( {x,t} \right)} = {\frac{\mathbb{d}\quad}{\mathbb{d}t}{{T\left( {x,t} \right)}.}}$Equation (17) is a conservation equation and is thus valid over theentire body, a part or a point of this body. Consequently, the integralsigns can be taken off: $\begin{matrix}{\underset{\underset{\underset{variation}{{Thermal}\quad{energy}}}{︸}}{c\quad{\rho(x)}\quad{h\left( {x,t} \right)}} = {\underset{\underset{\underset{entering}{{Rate}\quad{of}\quad{heat}}}{︸}}{{- \nabla} \cdot {q\left( {x,t} \right)}} + \underset{\underset{\underset{\underset{generation}{heat}}{{Rate}\quad{of}}}{︸}}{{\sigma\left( {x,t} \right)}.}}} & (18)\end{matrix}$

This equation Is said to be local, whereas Equation (17) is said to beglobal. The thermal energy variation is called the unsteady term, therate of heat entering is called the diffusion term and the rate of heatgeneration is called the source term.

Based on Equation (18), two cases are be considered: the steady statecase and the unsteady case with no source. The term steady simply meansthat there is no variation of the thermal energy of the system overtime. That is, the left side of Equation (18) is null:cρ(x)h(x, t)=0   (19)

This implies that the heat diffusion compensates for the internal heatproduction:∇·q(x, t)=σ(x, t)   (20)

In the unsteady case with no source we have:σ(x, t)=0   (21)which means that the time variation of thermal energy is explained bythe heat diffusion alone:cρ(x)h(x, t)=−∇·q(x, t)   (22)Constitutive Principles

In Equation (18), there are three unknown variables: ρ(x), h(x,t) andq(x,t). Let's look at the example of q(x,t). Suppose that we can measurethe time variation of the thermal energy (left side of the equation) andalso of the temperature T. We know that q(x,t) is related to thetemperature, but since different materials usually have differentdiffusion properties, the missing equation q(x,t)=f(T,x,t) must dependon properties of the material we are studying, such as its homogeneityand type of diffusivity. Consequently, the system of equations containsmore unknown variables than equations and the function f(T,x,t) must beadded to the system formed by Equation (18). This is due to thedifference in material properties. Different materials behavedifferently, but are subject to the same conservation laws. Constitutiveequations such as f(T,x,t), which reflect the internal constitution ofmaterials, allow us to complete the system of equations.

Decomposition Into Basic Laws

As indicated in the foregoing description, conservation equations arealways valid regardless of the materials, whereas constitutive equationsare dependent on their properties. When solving directly PDEs likeEquation (18) in a discretized context with methods such as the finitedifferences approach, one makes global assumptions about the time andspace behavior of the diffusion, energy variation and source termswithout taking into account the nature of the basic laws underlying theproblem. Some of these do not require any approximation since they comefrom conservation principles. Also, a more physically realistic solutioncan be obtained by choosing a proper approximation for each basic lawarising from a constitutive principle. Consequently, we propose todecompose the terms of Equation (18) into basic laws. This equation canbe broken down into one constitutive and two conservative laws for thesteady state case. In the unsteady case, an additional constitutive andanother conservative law must also be considered. Note that since thesource term is often known, it will not be decomposed. Recalling thatthe diffusion term α(x,t) is the rate of heat entering the particlelocated at x at time t, then:α(x, t)=−∇·q(x, t)   (23)is a first basic conservation law.

The second conservation law concerns the thermal tension. We firstdefine the thermal tension vector g(x,t) as the vector representing thedirection and magnitude of the greatest temperature decrease at a fixedtime t. As g(x,t) is source-oriented (from hot to cold), a minus signmust be inserted before ∇T(x,t) which represents the direction andmagnitude of the greatest temperature increase:g(x, t)=−∇T(x, t)   (24)This equation is a second basic law. Since the thermal tension is thegradient of a scalar field, it is by definition a conservative field inspace. It can be said that −T(x,t) is the potential field of g(x,t) [26,41].

The third law is a constitutive law. The heat flow density q(x,t) isdefined as the quantity and the direction of the heat flowing into theparticle located at point x at time t. It is represented by a vector andgreatly depends on the behavior of the material. In the case of uniform,homogeneous materials, it has been proven experimentally by Fourier [20,27] that q(x,t) is directly proportional to the difference oftemperature relative to neighbors of this particle:q(x, t)=λg(x, t)   (25)where λ is a material-specific thermal conductivity constant. The valueof λ is known for many types of materials. Equation (25) is called theFourier heat conduction law. For a non-homogeneous material, we considerthat it has the behavior of a homogeneous material on an infinitesimalpatch, but the conductivity changes with each patch; that is:q(x, t)=λ(x, t)g(x, t)

For the unsteady case, the fourth basic law is:ε(x, t)=cρ(x)h(x, t)   (26)where ε(x,t) is the thermal energy variation (the unsteady term). Thisequation is a constitutive one because it involves ρ(x), which ismaterial-dependent.

Finally, the fifth basic law is related to the temperature variation andis a conservation law: $\begin{matrix}{{h\left( {x,t} \right)} = {\frac{\mathbb{d}\quad}{\mathbb{d}t}{T\left( {x,t} \right)}}} & (27)\end{matrix}$

Considering only the temperature T_(x)(t) of the particle located at x,we reduce the basic law to a 1-dimensional equation and we can thus saythat T_(x)(t) is a conservative field in time-space.

To summarize, three basic laws for the diffusion term of equation (18)have been defined, that is:α(x, t)=−∇·q(x, t)q(x, t)=λg(x, t)g(x, t)=−∇T(x, t)

There are also two additional basic laws for the unsteady term, that is:ε(x, t)=cρ(x)h(x, t)${h\left( {x,t} \right)} = {\frac{\mathbb{d}\quad}{\mathbb{d}t}{T\left( {x,t} \right)}}$

Combining all these elements, the following relation is obtained:$\begin{matrix}{{c\quad{\rho\left( {x,t} \right)}\quad\frac{\mathbb{d}\quad}{\mathbb{d}t}{T\left( {x,t} \right)}} = {{\nabla{\cdot \left( {\lambda{\nabla{T\left( {x,t} \right)}}} \right)}} + {\sigma\left( {x,t} \right)}}} & (28)\end{matrix}$Discrete Representation of Images

Some algebraic tools used to model images will now be recalled from theabove description. An image is composed of two distinctive parts: theimage support (pixels) and some field quantity associated with eachpixel. This quantity may be scalar (e.g. gray level), vectorial (e.g.color, multispectral, optical flow) or tensorial (e.g. Hessian). Theimage support is modelled in terms of cubical complexes, chains andboundaries as described in the foregoing description. With theseconcepts, it is possible to give a formal description of an imagesupport of any dimension. For quantities, the concept of cochains hasbeen introduced, these cochains being representations of fields over acubical complex. For the use of these concepts in image processing, see[16].

As discussed hereinabove, an image is a complex of unit cubes usuallycalled pixels. A pixel γ ⊂ R^(n) is a product:γ=I ₁ ×I ₂ × . . . ×I _(n)where I_(j) is either a singleton or an interval of unit length withinteger endpoints. Thus I_(j) is either the singleton {k} and is said tobe a degenerate interval, or the closed interval [k, k+1] for some k εZ. The number q ε{0,1, . . . ,n} of non-degenerate intervals is bydefinition the dimension of γ, which is called a q-pixel. FIGS. 6 a-6 cillustrate three elementary pixels in R². For q≧1, let J={k₀,k₁, . . .,k_(q−1)} be the ordered subset {1,2, . . . ,n} of indices for whichI_(j) _(k) =[a_(j),b_(j)] is non-degenerate. Let us define:A _(k) _(j) σ=I ₁ × . . . I _(k) _(j) ⁻¹ ×{a _(j) }×I _(k) _(j) ₊₁ × . .. ×I _(n)andB _(k) _(j) σ=I ₁ × . . . I _(k) _(j) ⁻¹ ×{b _(j) }×I _(k) _(j) ₊₁ × . .. ×I _(n)

The A_(k) _(j) and the B_(k) _(j) are called the (q−1)-faces of σ. Onecan define the (q−2)-faces, . . . , down to the 0-faces of σ in the sameway. The faces of γ different from γ itself are called its proper faces.

By definition, a natural orientation of the cube is assumed for eachpixel. Suppose that γ denotes a particular positively oriented q-pixel.It is natural to denote the same pixel with opposite orientation by −γ.Examples of orientations are given in FIGS. 6 a-6 b. A cubical complexin R^(n) is a finite collection K of q-pixels such that every face ofany pixel of the image support is also a pixel in K and the intersectionof any two pixels of K is either empty or a face of each of them. Forexample, traditional 2D image models only considered pixels as 2D squareelements. The definitions presented above allow us to consider 2-pixels(square elements), 1-pixels (line elements) and 0-pixels (pointelements) simultaneously.

In order to write the image support in algebraic form, the concept ofchains is introduced. Any set of oriented q-pixels of a cubical complexcan be written in algebraic form by attributing to them the coefficient0,1 or −1, if they are not in the set or if they should or should not betaken with positive orientation, respectively. In order to representweighted domains, arbitrary integer multiplicity is allowed for eachq-pixel.

Given a topological space X ⊂ R^(n) in terms of a cubical complex, weget a free abelian group C_(q)(X) generated by all the q-pixels. Theelements of this group are called q-chains and they are formal linearcombinations of q-pixels [16]. A formal expression for a q-chain c_(q)is $c_{q} = {\sum\limits_{\gamma_{i} \in K}{\lambda_{i}\gamma_{i}}}$where λ_(i) ε Z.

The last step needed for the description of the image plane is theintroduction of the concept of a boundary of a chain. Given a q-pixel γ,we define its boundary ∂γ as the (q−1)-chain corresponding to thealternating sum of its (q−1)-faces. The sum is taken according to theorientation of the (q−1)-faces with respect to the orientation of theq-pixel. A (q−1)-face of γ is said to be positively oriented relative tothe orientation of γ if its orientation is compatible with theorientation of γ. By linearity, the extension of the definition ofboundary to arbitrary q-chains is easy. For instance, in FIGS. 6 b and 6c, the boundary of the 1-pixel a is x₂−x₁ and the boundary of the2-pixel A is a+b−c−d; then a and b are said to be positively orientedwith respect to the orientation of A but c and d are said to benegatively oriented with respect to the orientation of A. Let us noticethat the boundary of a 1-pixel is always the difference between itsboundary points. The boundary can be defined recursively. Given a(q−1chain and a q-chain γ_(q) defined as γ_(q)=γ_(q−1)×[a,b], theboundary of γ_(q) can be recursively written as:∂γ_(q)=∂γ_(q−1) ×[a, b]+(−1)^((q−1))(γ_(q−1) ×{b}−γ _(q−1) ×{a})   (29)

In order to model the pixel quantity over the image plane, anapplication F that associates a global quantity with all q-pixels γ of acubical complex is determined and is denoted by <F,γ>. This quantity maybe any mathematical entity such as a scalar, a vector, etc. For twoadjacent q-pixels γ₁ and γ₂, F must satisfy <F, λ₁γ₁+λ₂γ₂>=λ₁<F,γ₁>+λ₂<F, γ₂>, which means that the sum of the quantity over each pixelis equal to the quantity over the two pixels. The resultingtransformation F: C_(q)(X)→R is called a q-cochain and is used as arepresentation of a quantity over the cubical complex.

Finally, an operator is needed to associate a global quantity with the(q+1)-pixels according to the global quantities given on their q-faces.Given a q-cochain F, we define an operator ∂, called the coboundaryoperator, which transforms F into a (q+1)-cochain ∂F such that:<δF, γ>=<F, ∂γ>  (30)for all (q+1 )-chains γ. The coboundary is defined as the signed sum ofthe physical quantities associated with the q-faces of γ. The sum istaken according to the relative orientation of the q-faces of the (q+1)-pixels of γ with respect to the orientation of the pixels. FIG. 7presents an example of the coboundary operator for a 2-pixel.

With this image model in hand, the basic laws described hereinabove willbe used to rewrite the global heat transfer equation in algebraic terms[43].

Representation of the Heat Transfer Equation

The process for representing the heat transfer equation in terms ofalgebraic topology can be summarized as follows. The image support issubdivided into cubical complexes. Basic laws are applied to pixels ofvarious dimensions. These laws involve the computation of globalquantities on pixels, expressed as cochains. Some of these laws linkglobal quantities on pixels with, global quantities on their boundariesand hence are expressed as coboundaries. The other laws are expressed aslinear transformations between pairs of cochains. The topologicalformalism of cochain and coboundary is a generic one; that is, it doesnot offer computational rules. The cochains must be instantiateddepending on the problem to be considered.

The basic laws presented hereinabove will be reformulated in atopological way and then to give computational rules for cochains in thecontext of the heat transfer problem. Since we want to represent twokinds of global values over the spatio-temporal image, two complexeswill be used. The first complex is associated with global valuescorresponding to projections onto the tangential part of the domain(e.g. global thermal tension) while the second complex refers to valuesrelated to projections onto the normal part of the domain (e.g. heatentering a particle). These two distinct orientations (see FIGS. 8 a and8 b give rise to two different complexes.

Global Heat Transfer

Let us assume that an image has n spatial dimensions and r pixels.Suppose also that a time interval [t₀, t₁] can be split into i equalsub-intervals [t_(k), t_(k+1)], k ε [0, i−1]. Let us consider ann-complex representing the subdivided spatial support of the imageK^(s′). One can consider an (n+1)-complex representing thespatio-temporal support of the image:K ^(s) =K ^(sl) ×[t _(k) , t _(k+1) ], ∀k ε [0, l−1]

Now, let us consider γ_(E), an (n+1)-pixel of K^(s), and the followingcochain, defined as <ε,γ_(E)>. Thus, we therefore need to define whichvalue to use as cochain ε in the heat transfer problem. Let us define εas the global energy variation of the (n+1)-pixel γ_(E): $\begin{matrix}{\left\langle {\mathcal{E},\gamma_{E}} \right\rangle = {\int_{\gamma_{E}}^{\quad}{{\varepsilon\left( {x,t} \right)}\quad{\mathbb{d}\gamma_{E}}}}} & (31)\end{matrix}$where dγ_(E) is an infinitesimal part of the domain represented byγ_(E). Now, using the global version of Equation (18), we obtain:∫_(γ_(E))  ε(x, t)  𝕕γ_(E) = ∫_(γ_(E))  α(x, t)  𝕕γ_(E) + ∫_(γ_(E))  σ(x, t)  𝕕γ_(E)From this equation, we define two more cochains, representing first, theglobal diffusion: $\begin{matrix}{\left\langle {\mathcal{D},\gamma_{E}} \right\rangle = {\int_{\gamma_{E}}^{\quad}{{\alpha\left( {x,t} \right)}\quad{\mathbb{d}\gamma_{E}}}}} & (32)\end{matrix}$and second, the global source: ⟨𝒮, γ_(E)⟩ = ∫_(γ_(E))  σ(x, t)  𝕕γ_(E)

Thus, the following relation is obtained between the three cochains:<ε, γ_(E) >=<D, γ _(E) >+<S, γ _(E)>  (33)

The rules used for cochains ε and D are then decomposed into basic laws.The rule for cochain S is not decomposed since it is assumed that itsglobal value is known on γ_(E). Let us finally mention that both steadyand unsteady heat transfer problems can be considered using Equation(33) by setting respectively, cochains ε and S, to zero.

Global Temperature Variation

Let us consider another n-complex, K^(p′), representing the subdividedspatial domain of the image. An (n+1 )-complex representing thespatio-temporal image can then be defined as:K ^(p) =K ^(p′) ×[t _(k) , t _(k+1) ], ∀k ε [0, l−1]

Let us consider γ_(H), a 1-pixel of K^(p) defined asx_(i)×[t_(k),t_(k+1)], i ε[1,r], kε[0,l−1] where x_(i) is a 0-pixel ofK^(p′). Let us also consider a 0-cochain T and a 1-cochain H such that:<H, γ_(H)>=<δT, γ_(H)>=<T, ∂γ_(H)>  (34)

FIGS. 9 a and 9 b present examples of cochains T and H for K^(p) ofdimension 3.

Applying Equation (29), it is found that the boundary of γ_(H) is∂γ_(H)=x_(i)×{t_(k+1)}−x_(i)×{t_(k)}. According to the linearity of thecochain, the computational rule relating the global value associated toγ_(H) with the values at its boundary x_(i)×{t_(k)} and x_(i)×{t_(k+1)}is:<T, ∂γ _(H) >=<T, x _(i) ×{t _(k+1) }−x _(i) ×{t _(k) }>=<T, x _(i) ×{t_(k+1) }>−<T, x _(i) ×{t _(k)}>  (35)

This equation is general and applies to many problems. To define whichvalues to use as 0-cochain and 1-cochain, let us take the global versionof Equation (27) on γ_(H) and apply the fundamental calculus theorem:${\int_{\gamma_{H}}^{\quad}{{h\left( {x,t} \right)}\quad{\mathbb{d}\gamma_{H}}}} = {{\int_{t_{k}}^{t_{k + 1}}{\frac{\mathbb{d}\quad}{\mathbb{d}t}{T\left( {x_{i},t} \right)}\quad{\mathbb{d}t}}} = {{T\left( {x_{i},t_{k + 1}} \right)} - {T\left( {x_{i},t_{k}} \right)}}}$

Looking at this equation, it can be seen that it is similar to Equation(35). Thus we define T=T(x, t). The location of the unknown temperaturesto compute will correspond to the 0-pixels of K^(p). In order to fulfillEquation (34), the following relation is used: $\begin{matrix}{\left\langle {\mathcal{H},\gamma_{H}} \right\rangle = {\int_{\gamma_{H}}^{\quad}{{h\left( {x,t} \right)}\quad{\mathbb{d}\gamma_{H}}}}} & (36)\end{matrix}$this relation being called the global temperature variation. These threeequations are extended by linearity to a 1-chain of K^(p) defined asγ×[t_(k), t_(k+1)], where γ is an arbitrary 0-chain of K^(p′).Global Energy Variation

We want to link cochains H and ε, representing the global temperaturevariation and the global energy variation, respectively. For thispurpose, a representation of Equation (26) is needed. The two cochainsare not from the same cubical complex (H is from K^(p) and ε is fromK^(s)), and moreover, Equation (26) is material-dependent; thereforethey cannot be linked exactly. However, we can express this link as alinear transformation:H→ε

Recalling Equation (31), the following relation is obtained:$\begin{matrix}{\left\langle {\mathcal{E},\gamma_{E}} \right\rangle = {{\int_{\gamma_{E}}^{\quad}{{\varepsilon\left( {x,t} \right)}\quad{\mathbb{d}\gamma_{E}}}} = {\int_{\gamma_{E}}^{\quad}{c\quad{\rho(x)}\quad{h\left( {x,t} \right)}\quad{\mathbb{d}\gamma_{E}}}}}} & (37)\end{matrix}$

Unfortunately, the value of ρ(x) or h(x, t) is not known at all pointsof the volume. Consequently, these two fields are approximated over thevolume. The approximations are denoted by {tilde over (ρ)}(x) and {tildeover (h)}(x, t). For one 1-pixel γ_(H), defined asx_(i)×[t_(k),t_(k+1)], the approximation is performed piecewise suchthat {tilde over (h)}(x, t) must fulfill Equation (36). $\begin{matrix}{{\int_{t_{k}}^{t_{k + 1}}{\overset{\sim}{h}\left( {x_{i},t} \right)}} = \left\langle {\mathcal{H},\gamma_{H}} \right\rangle} & (38)\end{matrix}$

Equation (37) thus becomes: $\begin{matrix}{\left\langle {\mathcal{E},\gamma_{E}} \right\rangle = {{\int_{\gamma_{E}}^{\quad}{c\quad{\overset{\sim}{\rho}(x)}\quad{\overset{\sim}{h}\left( {x,t} \right)}\quad{\mathbb{d}\gamma_{E}}}} = {{f_{e}\left( {c,\mathcal{H}} \right)} = \Gamma}}} & (39)\end{matrix}$where dV is an infinitesimal part of γ_(E) which depends on the choiceof the approximation functions {tilde over (ρ)}(x) and {tilde over(h)}(x, t) and the position of K^(s) with respect to K^(p).Global Diffusion

Let us consider an n-cochain Q and an (n+1)-cochain D defined by thecoboundary:<D, γ_(E)>=<δQ, γ_(E)>=<Q, ∂γ_(E)>  (40)

FIG. 10 presents examples of Q and D for K^(s) of dimension 3. Let usassume that the n-faces γ_(Q) _(i) of γ_(E) are positively orientedrelative to γ_(E). According to the linearity of the cochain, thecomputational rule is: $\begin{matrix}{\left\langle {\mathcal{D},\gamma_{E}} \right\rangle = {\sum\limits_{\gamma_{Q_{i}} \in {\partial\gamma_{E}}}\left\langle {Q,\gamma_{Q_{i}}} \right\rangle}} & (41)\end{matrix}$

Again, this equation is general; hence a global value is found for the(n+1 )-cochain D, which can be computed by summing the global values atthe boundary of γ_(E). According to Equation (32), the followingrelation is obtained: ⟨𝒟, γ_(E)⟩ = ∫_(γ_(E))  α(x, t)  𝕕γ_(E)

The divergence theorem is applied to this equation to obtain:${{\mathcal{g}}_{t}(x)} = {\frac{1}{4\quad\pi\quad t}{\mathbb{e}}^{- {(\frac{x^{2} + y^{2}}{4t})}}}$where n(x, t) is the normal vector to an infinitesimal part of thedomain represented by γ_(Q). This last equation is in the form of acoboundary (Equation (41)), from the following relation is defined:$\begin{matrix}{\left\langle {Q,\gamma_{Q_{i}}} \right\rangle = {\int_{\gamma_{Q_{i}}}^{\quad}{{{- {q\left( {x,t} \right)}} \cdot {n\left( {x,t} \right)}}\quad{\mathbb{d}\gamma_{Q_{i}}}}}} & (42)\end{matrix}$

Again, the previous definitions can be extended by linearity toarbitrary (n+1 )-chains of K^(s). And there is absolutely noapproximation in these equations.

Global Thermal Tension

Let us consider a 1-pixel γ_(G) of K^(p) defined as γ×t_(k), k ε80,l−1], where γ is a 1-pixel of K^(p′) whose boundary is defined as∂γ=x_(j)−x_(i), i,j ε[1,r]. Let us also consider a 0-cochain T and a1-cochain G defined by the coboundary:<G, γ_(G)>=<δT, γ_(G)>=<T, ∂γ_(G)>  (43)

FIGS. 11 a and 11 b present examples of cochains T and G for one 3-pixelof K^(p). The boundary of γ_(G) is ∂γ_(G)=x_(j)×{t_(k)}−x_(i)×{t_(k)}.According to the linearity of the cochain, the computational rulerelating the global value associated with γ_(G) to the values atx_(i)×{t_(k)} and x_(j)×{t_(k)} is:<T, ∂γ _(G) >=<T, x _(j) ×{t _(k) }−x _(i) ×{t _(k) >=<T, x _(j) ×{t_(k) }>−<T, x _(i) ×{t _(k)}>  (44)

To define which values to use as cochains G and T let us take the globalform of Equation (24) on γ_(G): $\begin{matrix}\begin{matrix}{{\int_{\gamma_{G}}^{\quad}{{g\left( {x,t} \right)} \cdot {\mathbb{d}\gamma_{G}}}} = {\int_{x_{i}}^{x_{j}}{{g\left( {x,t_{k}} \right)} \cdot {\mathbb{d}\gamma}}}} \\{\quad{= {\int_{x_{i}}^{x_{j}}{{- {\nabla{T\left( {x,t_{k}} \right)}}} \cdot {\mathbb{d}\gamma}}}}}\end{matrix} & (45)\end{matrix}$where dγ is an infinitesimal part of γ. Since g(x, t) is a spatialconservative field, we can apply the line integral theorem [26, 41]saying that for a conservative field F(x)=∇f(x) and two points A and B,in an open connected region containing F(x), the integral of thetangential part of F(x) along the curve R joining A and B is independentof the path (FIG. 12): $\begin{matrix}{{\int_{A}^{B}{{F(x)} \cdot {\mathbb{d}R_{1}}}} = {\int_{A}^{B}{{F(x)} \cdot {\mathbb{d}R_{2}}}}} \\{\quad{= {\int_{A}^{B}{{F(x)} \cdot {\mathbb{d}R_{3}}}}}} \\{\quad{= {{f(B)} - {f(A)}}}}\end{matrix}$

From this theorem, Equation (45) can be rewritten as: $\begin{matrix}\begin{matrix}{{\int_{x_{i}}^{x_{j}}{{g\left( {x,t_{k}} \right)} \cdot {\mathbb{d}\gamma}}} = {\left( {- {T\left( {x_{j},t_{k}} \right)}} \right) - \left( {- {T\left( {x_{i},t_{k}} \right)}} \right)}} \\{\quad{= {{T\left( {x_{i},t_{k}} \right)} - {T\left( {x_{j},t_{k}} \right)}}}}\end{matrix} & (46)\end{matrix}$

Looking at Equation (46), it can be seen that it is similar to Equation(44). Thus T=T(x, t) is defined. Consequently, the location of theunknown temperatures to be computed correspond to the 0-pixels of K^(p)which is coherent with the conclusions hereinabove. In order to fulfillEquation (43), we have: $\begin{matrix}{\left\langle {\mathcal{G},\gamma_{G}} \right\rangle = {\int_{\gamma_{G}}^{\quad}{{- {g\left( {x,t} \right)}} \cdot {\mathbb{d}\gamma_{G}}}}} & (47)\end{matrix}$

The previous definitions are extended by linearity to 1-chains of K^(p)defined as γ×{t_(k)}, where γ is an arbitrary 1-chain of K^(p′).

Heat Flow Density

The coboundaries <Q, ∂γ_(E)> (Equation (40)) and <T, ∂γ_(G)> (Equation(43)) provide exact global versions of Equation (23) on K^(s) andEquation (24) on K^(p), respectively. In order to complete the diffusionterm, Equation 25, which links local values g(x,t) and q(x,t) isrepresented. Equation (25) is a constitutive equation and cannot berepresented by a topological equation. However, a relation transformingcochain G into cochain Q can be found:<G, γ_(G)>

<Q, γ_(Q)>as a global counterpart for Equation (25). To find this transformation,Equation 42 is recalled: $\begin{matrix}{\left\langle {Q,\gamma_{Q_{i}}} \right\rangle = {\int_{\gamma_{Q_{i}}}^{\quad}{{{- {q\left( {x,t} \right)}} \cdot {n\left( {x,t} \right)}}\quad{\mathbb{d}\gamma_{Q_{i}}}}}} \\{\quad{= {\int_{\gamma_{Q_{i}}}^{\quad}{{- \lambda}\quad{{g\left( {x,t} \right)} \cdot {n\left( {x,t} \right)}}\quad{\mathbb{d}\gamma_{Q_{i}}}}}}}\end{matrix}$this equation relating cochain Q to field g(x,t). Unfortunately, fieldg(x,t) is not known, so that this equation has to be approximated with afield {tilde over (g)}(x,t). Let us consider γ_(n), an n-pixel of K^(p)defined as γ_(n)=γ_(x)×{t_(k)}, k ε [O,l−1] where γ_(x) is an n-pixel ofK^(p′). This approximation is performed piecewise such that for each1-face γ_(G) of y_(n), {tilde over (g)}(x,t) satisfies: $\begin{matrix}{{\int_{\gamma_{G}}^{\quad}{{- {\overset{\sim}{g}\left( {x,t} \right)}} \cdot {\mathbb{d}R}}} = \left\langle {\mathcal{G},\gamma_{G}} \right\rangle} & (48)\end{matrix}$where dR is an infinitesimal part of the domain represented by γ_(G).Equation (25) is then applied to obtain {tilde over (q)}(x,t):{tilde over (q)}(x, t)=λ{tilde over (g)}(x, t)at all points of the domain. Equation (42) becomes: $\begin{matrix}{\left\langle {Q,\gamma_{Q_{i}}} \right\rangle = {{\int_{\gamma_{Q_{i}}}^{\quad}{{{- {\overset{\sim}{q}\left( {x,t} \right)}} \cdot {n\left( {x,t} \right)}}\quad{\mathbb{d}\gamma_{Q_{i}}}}} = {f_{\mathcal{g}}\left( {\lambda,\mathcal{G}} \right)}}} & (49)\end{matrix}$

The transformation that is looked for is thus:Λ=f _(g)(λ, G)which depends on the choice of an approximation function {tilde over(g)}(x,t) and the position of K^(s) with respect to K^(p).Boundary Conditions

The decomposition process that has been presented herein is carried outwith the assumption that all the needed quantities surrounding a pixelare known. For instance, in the-steady state heat transfer problem, fora particular (n+1)-pixel, the cochain S is known for all othersurrounding (n+1)-pixels, that is, there are as many equations asvariables.

Unfortunately, this assumption is not verified at the borders of theimage. Thus, as in solving the PDE, certain boundary conditions areimposed to specify the gray-level conditions at the boundary of theimage. For instance, these conditions may prescribe the values of eithercochain T (Dirichlet boundary conditions) or cochain Q (Neumann boundaryconditions).

Summary of the Algorithm

The algorithm used to find an expression of the temperatures at timet_(k+1) as a function of the temperatures at time t_(k) will now besummarized. The input data for this algorithm are the cochain S and theDirichlet boundary conditions. That is, T is known for all pixels on theboundary of the image, which includes the values at time t₀.

-   -   1. Choose the positions for K^(p′) and K^(s′).    -   2. Compute ε as a function of H:        -   (a) Choose the approximation functions {tilde over (h)}(x,t)            and {tilde over (ρ)}(x).        -   (b) Apply Equation (38) to find {tilde over (h)}(x,t_(k),            t_(l)) as a function of H.        -   (c) Apply Equation (39) to find the transformation Γ,            expressing ε as a function of H.    -   3. Apply Γ and Equation (34) to find ε as a function of T.    -   4. Compute Q as a function of G:        -   (a) Choose the approximation function {tilde over (g)}(x,t).        -   (b) Apply Equation (48) to find {tilde over (g)}(x,t) as a            function of G.        -   (c) Apply Equation (49) to find the transformation Λ,            expressing Q as a function of G.    -   5. Apply Equation (40), Λ and Equation (43) to find D as a        function of T.    -   6. Apply Equation (33) to obtain an equation of the temperatures        at time t_(k+1) as a function of the temperatures at time t_(k)

FIG. 13 presents an overview of the computational scheme.

Applications

Linear Isotropic Diffusion

One of the most direct applications of the heat transfer equation is theisotropic diffusion of gray-level intensities; that is, smoothing. For a2D image I(x), with x=(x,y), the resolution of the PDE: $\begin{matrix}{{\frac{\partial}{\partial t}{I\left( {x,t} \right)}} = {\nabla^{2}{I\left( {x,t} \right)}}} & (50)\end{matrix}$is equivalent to the convolution:I(x, t)=(I*g _(t))(x)where${{\mathcal{g}}_{t}(x)} = {\frac{1}{4\quad\pi\quad t}\quad{\mathbb{e}}^{- {(\frac{x^{2} + y^{2}}{4t})}}}$is a Gaussian with variance σ²=2t [25]. One can see t as the scale ofthe smoothing operation. Let us assume that the Laplacian image at scalet:L(x, t)=∇² I(x, t)   (51)is known. One can consider this equation as a steady state heat transferproblem with T(x,t)=I(x,t), σ(x,t)=−L(x,t) and λ=1.

It is desired to solve Equation (51) for local I(x,t) located at thecenter of each image pixel. Employing the process presented hereinabove,we first position the two cubical complexes representing twosubdivisions of the image plane. As stated hereinabove, the primarycomplex K^(p′) is defined with 0-pixels corresponding to pixel centers.For the sake of simplicity, K^(s′) corresponds to the image pixels; thatis, the secondary 2-pixels γ_(s) are rectangular and symmetricallystaggered relative to the 1-pixels of K^(p) and the 1-pixels γ_(q) ofK^(s) intersect orthogonally in the centers of the primary 1-pixels.Since there is no variation in steady-state heat transfer over time, thetime parameter is dropped. This means that K^(p)=K^(p′), K^(s)=K^(s′)and the time integral in cochain computation are dropped. It can be seenthat the approximation function f_(g) depends on the position of K^(s)with respect to K^(p).

FIG. 14 shows the two complexes for a 5×5 image. Positioning the 1-facesof K^(s) such as each passes through the center point between two0-pixels of K^(p) allows us to compute a polynomial function of order 1with the same accuracy as that obtained using one of order 2 [37, 34].

A global value for the 2-cochain <S, γ_(E)> is needed. If it is assumedthat a pixel value represents the global value of intensity, <S,γ_(E)>=−L(x) can be directly set. This assumption is reasonable if imageacquisition is considered as a process which accumulates the totalnumber of photons within a global area corresponding to the pixel [22].

An approximation function {tilde over (g)}(x) is chosen. For simplicity,we assume that {tilde over (g)}(x) arises from a bilinear approximation,that is:{tilde over (g)}(x)=(a+by)·{right arrow over (i)}+(c+dx)·{right arrowover (j)}

Given a 2-pixel γ_(p) of K^(p), {tilde over (g)}(x) satisfies Equation48 for each 1-face of γ_(p). As an example, let us find the coefficientsa, b, c and d for such a pixel defined as in FIG. 15: $\begin{matrix}{\mathcal{G}_{1} = {\int_{0}^{\Delta}{{{- {\overset{\sim}{g}\left( {x,0} \right)}} \cdot \overset{\rightarrow}{i}}\quad{\mathbb{d}x}}}} \\{\mathcal{G}_{2} = {\int_{0}^{\Delta}{{{- {\overset{\sim}{g}\left( {\Delta,y} \right)}} \cdot \overset{\rightarrow}{j}}\quad{\mathbb{d}y}}}} \\{\mathcal{G}_{3} = {\int_{0}^{\Delta}{{{- {\overset{\sim}{g}\left( {x,\Delta} \right)}} \cdot \overset{\rightarrow}{i}}\quad{\mathbb{d}x}}}} \\{\mathcal{G}_{4} = {\int_{0}^{\Delta}{{{- {\overset{\sim}{g}\left( {0,y} \right)}} \cdot \overset{\rightarrow}{j}}\quad{\mathbb{d}y}}}}\end{matrix}$from which $\begin{matrix}\begin{matrix}{{{\overset{\sim}{g}(x)} = {- {\frac{1}{\Delta}\left\lbrack {{\left( {\mathcal{G}_{1} + {\frac{\left( {\mathcal{G}_{3} - \mathcal{G}_{1}} \right)}{\Delta}y}} \right) \cdot \overset{\rightarrow}{i}} + {\left( {\mathcal{G}_{4} + {\frac{\left( {\mathcal{G}_{2} - \mathcal{G}_{4}} \right)}{\Delta}x}} \right) \cdot \overset{\rightarrow}{j}}} \right\rbrack}}},} \\{x \in \gamma_{p}}\end{matrix} & (52)\end{matrix}$is obtained. {tilde over (g)}(x) is thus a piecewise function of G, butas G is computed from T, and {tilde over (g)}(x) can also be expressedas a function of T. For each primary 2-pixel, Equation 25 is applied toobtain {tilde over (q)}(x)={tilde over (g)}(x).

The next step is to compute <Q, γ_(Q)> for K^(s) from Equation (49).Each secondary 2-pixel γ_(s) intersects with four primary 2-pixels,γ_(pa), γ_(pb), γ_(pc) and γ_(pd). There are four segments in theapproximation function {tilde over (q)}(x) corresponding to the fourprimary 2-pixels; that is, {tilde over (q)}_(a)(x), {tilde over(q)}_(b)(x), {tilde over (q)}_(c) (x) and {tilde over (q)}_(d)(x). FIGS.16 a-16 c illustrate γ_(E). Cochain <Q, γ_(Q)> corresponding to the four1-faces of γ_(E) is found by: $\begin{matrix}\begin{matrix}{Q_{1} = {{\int_{0}^{\Delta/2}{{{- {{\overset{\sim}{q}}_{a}\left( {{\Delta/2},y} \right)}} \cdot \overset{\rightarrow}{i}}\quad{\mathbb{d}y}}} + {\int_{{- \Delta}/2}^{0}{{{- {{\overset{\sim}{q}}_{b}\left( {{\Delta/2},y} \right)}} \cdot \overset{\rightarrow}{i}}\quad{\mathbb{d}y}}}}} \\{\quad{= {\frac{3\quad\mathcal{G}_{1}}{4} + \frac{\mathcal{G}_{3}}{8} + \frac{\mathcal{G}_{5}}{8}}}} \\{Q_{2} = {{\int_{0}^{\Delta/2}{{{- {{\overset{\sim}{q}}_{b}\left( {x,{{- \Delta}/2}} \right)}} \cdot \left( {- \overset{\rightarrow}{j}} \right)}\quad{\mathbb{d}x}}} +}} \\{\quad{\int_{{- \Delta}/2}^{0}{{{- {{\overset{\sim}{q}}_{c}\left( {x,{{- \Delta}/2}} \right)}} \cdot \left( {- \overset{\rightarrow}{j}} \right)}\quad{\mathbb{d}x}}}} \\{\quad{= {{- \frac{3\quad\mathcal{G}_{7}}{4}} - \frac{\mathcal{G}_{6}}{8} - \frac{\mathcal{G}_{9}}{8}}}} \\{Q_{3} = {{\int_{0}^{\Delta/2}{{{- {{\overset{\sim}{q}}_{a}\left( {x,{\Delta/2}} \right)}} \cdot \overset{\rightarrow}{j}}\quad{\mathbb{d}x}}} + {\int_{{- \Delta}/2}^{0}{{{- {{\overset{\sim}{q}}_{d}\left( {x,{\Delta/2}} \right)}} \cdot \overset{\rightarrow}{j}}\quad{\mathbb{d}x}}}}} \\{\quad{= {\frac{3\quad\mathcal{G}_{4}}{4} + \frac{\mathcal{G}_{2}}{8} + \frac{\mathcal{G}_{11}}{8}}}} \\{Q_{4} = {{\int_{0}^{\Delta/2}{{{- {{\overset{\sim}{q}}_{d}\left( {{{- \Delta}/2},y} \right)}} \cdot \left( {- \overset{\rightarrow}{i}} \right)}\quad{\mathbb{d}y}}} +}} \\{\quad{\int_{{- \Delta}/2}^{0}{{{- {{\overset{\sim}{q}}_{c}\left( {{{- \Delta}/2},y} \right)}} \cdot \left( {- \overset{\rightarrow}{i}} \right)}\quad{\mathbb{d}y}}}} \\{\quad{= {{- \frac{3\quad\mathcal{G}_{10}}{4}} - \frac{\mathcal{G}_{12}}{8} - \frac{\mathcal{G}_{8}}{8}}}}\end{matrix} & (53)\end{matrix}$

Using Equation 41, we obtain:<D, γ _(E) >=Q ₁ +Q ₂ +Q ₃ +Q ₄   (54)

Substituting Equation (43) in Equation (53), Equation (53) in Equation(54) and Equation (54) in Equation (33), <S, γ_(E)> can now be expressedas a function of T. As an example, <S, γ_(E)> is presented for a 2-pixelγ_(E), and defined as in FIG. 16 a-16 c: $\begin{matrix}\begin{matrix}{\left\langle {\mathcal{S},\gamma_{s}} \right\rangle = {{{- 3}\quad\mathcal{T}_{0,0}} + {\frac{1}{2}\left\lbrack {\mathcal{T}_{0,1} + \mathcal{T}_{1,0} + \mathcal{T}_{0,{- 1}} + \mathcal{T}_{{- 1},0}} \right\rbrack} +}} \\{\frac{1}{4}\left\lbrack {\mathcal{T}_{{- 1},1} + \mathcal{T}_{1,1} + \mathcal{T}_{1,{- 1}} + \mathcal{T}_{{- 1},{- 1}}} \right\rbrack}\end{matrix} & (55)\end{matrix}$

For each non-border pixel (represented by a secondary 2-pixel), anequation in the form of Equation (55) is obtained. For the-borderpixels, T=I(x) is set. Solving this system, the smoothed image I(x,t)=Tis obtained.

Optical Flow

An indirect application of the heat transfer equation is the computationof optical flow for a 2D image sequence I(x,t), using the Horn andSchunk [29] algorithm. It can be shown that the velocity vectoru(x,t)=(u(x,t), v(x,t)) satisfies the following constraint arising fromvariational calculus (for greater legibility, (x,t) has been dropped):I _(x) ² u+I _(x) I _(y) v=α ²∇² u−I _(x) I _(t)I _(x) I _(y) u+I _(y) ² v=α ²∇² v−I _(y) I _(t)   (56)where α is a weighting factor and I_(x), I_(y) and I_(t) are the firstderivatives of I(x,t) in x, y and t, respectively. Let us rewriteEquation (56) in the following vectorial form:∇I(∇I·u)=α²∇² u−I _(t) ∇Iwhere ∇²u=(∇²u, ∇²v). Reorganizing the terms of Equation (56), thefollowing equation is obtained:α²∇² u=∇I(∇I·u)+I _(t) ∇I   (57)

Taking σ(x, t)=−∇I(∇I·u)−I_(t) ∇I as a heat source, Equation (57) can beseen as a steady state heat transfer equation in which the cochain Tcorresponds to u(x, t) and λ=α². It can thus be decomposed using themethod described hereinabove. The cochain T is u(x, t), and thefollowing relation is obtained: $\begin{matrix}{{- \frac{I_{t}{\nabla I}}{\alpha^{2}}} = {{{- 3}\quad\mathcal{U}_{0,0}} + {{\nabla I}\left( {{\nabla I} \cdot \mathcal{U}_{0,0}} \right)} +}} \\{{\frac{1}{2}\left\lbrack {\mathcal{U}_{0,1} + \mathcal{U}_{1,0} + \mathcal{U}_{0,{- 1}} + \mathcal{U}_{{- 1},0}} \right\rbrack} +} \\{\frac{1}{4}\left\lbrack {\mathcal{U}_{{- 1},1} + \mathcal{U}_{1,1} + \mathcal{U}_{1,{- 1}} + \mathcal{U}_{{- 1},{- 1}}} \right\rbrack}\end{matrix}$

For the same reasons as in the linear diffusion problem, specialconsiderations are needed at the borders of the image. Zero velocity isassumed at the borders of the image and the system is solved to get thevelocity field for each point of the image.

Nonlinear Diffusion

Linear isotropic diffusion reduces noise but also blurs edges. As thescale increases, edges tend to be harder to identify [43]. One possibleway of reducing this effect might be to consider the heat conductioncoefficient λ as a field function dependent on the magnitude of theedges; that is: $\begin{matrix}{{\frac{\partial}{\partial t}{I\left( {x,t} \right)}} = {\nabla{\cdot \left( {{{\mathcal{g}}\left( {{\nabla{I\left( {x,t} \right)}}}^{2} \right)}\quad{\nabla{I\left( {x,t} \right)}}} \right)}}} & (58)\end{matrix}$which corresponds to Equation (28) with λ(x,t)=g(|∇I(x,t)|²),T(x,t)=I(x,t), ρ(x)=1, c=1 and σ(x,t)=0 (i.e., unsteady transfer with nosource). The conduction function g(s) must display the followingbehavior: in constant regions, there should be linear isotropicdiffusion (Equation (50)), that is g(|∇I(x,t)|²)=1 for |∇I(x,t)|²=0, andalmost no diffusion when the magnitude of the edge is great; that is,g(|∇I(x,t)|²)=0 for |∇I(x,t)|²→∞. Perona and Malik [38] proposed thefollowing functions: $\begin{matrix}{{{{\mathcal{g}}(s)} = \frac{1}{1 + \frac{s^{2}}{k^{2}}}},} & \left( {k > 0} \right) \\{and} & \quad \\{{{{\mathcal{g}}(s)} = {\mathbb{e}}^{- \frac{s^{2}}{k^{2}}}},} & \left( {k > 0} \right)\end{matrix}$

The parameter k in these functions is difficult to set because itcontrols the threshold of diffusion but also the steepness of thefunction [35]. An advantageous alternative is to use the function:${{\mathcal{g}}(s)} = {\frac{1}{2}{\tanh\left( {{\gamma\left( {k - s} \right)} + 1} \right)}}$where k and γ control the threshold and the steepness, respectively.Equation (58) is then soved for a particular t (the scale) with initialconditions:I(x, 0)=I(x)where I(x) is the original image.

Let us assume that we have I time steps Δt=t/l. First, the same cubicalcomplexes K^(p′) and K^(s′) are used as hereinabove and the followingrelations are defined:K ^(p) =K ^(p′) ×[t _(k) , t _(k+1)]K ⁸ =K ^(8′) ×[t _(k) , t _(k+1) ], t _(k) =kΔt, ∀k ε [0, l−1]

Secondly, the following assumption is made about the spatial behavior ofh(x, t); that is, the approximation function {tilde over (h)}(x, t) ischosen. For a 3-pixel γ_(E) as defined in FIG. 17, it is assume that His the mean value over [−Δ/2, Δ/2]×[−Δ/2, Δ/2]. Thus, ε=H; that is,using Equation (34), the following relation can also be written:<ε, γ_(E) >=T ¹ −T ⁰   (59)

For the sake of simplicity, the same spatial bilinear approximationfunction {tilde over (g)}(x, t) as hereinabove is used. The behaviorover a time step has to be approximated. Some common assumptions abouttime variation may be generalized by proposing [37]: $\begin{matrix}{{{\int_{t_{k}}^{t_{k + 1}}{{A(t)}\quad{\mathbb{d}t}}} = {\left( {{w\quad{A\left( t_{k + 1} \right)}} + {\left( {1 - w} \right){A\left( t_{k} \right)}}} \right)\quad\Delta\quad t}},} & {0 \leq w \leq 1}\end{matrix}$where A(t) is some quantity and w is a weighting factor [37]. Somevalues of w lead to well-known schemes: 1) w leads to the explicitscheme; that is, the value at t_(k) prevails for the entire timeinterval except at time t_(k+1), 2) w=1 leads to the fully implicitscheme; that is, the value changes at time t_(k) from A(t_(k)) toA(t_(k+1)) and stays there throughout the whole time interval. 3) w=0.5leads to the semi-implicit or Crank-Nicolson scheme; that is, there is alinear variation of A(t). It is proposed to use the implicit schemebecause for large values of Δt, it best emulates long term time behaviorfor heat [37]. That is for a 3-pixel, and for w=1: $\begin{matrix}{{{\overset{\sim}{g}(x)} = {{- {\frac{1}{\Delta}\left\lbrack {{\left( {\mathcal{G}_{1}^{1} + {\frac{\left( {\mathcal{G}_{3}^{1} - \mathcal{G}_{1}^{1}} \right.}{\Delta}y}} \right) \cdot \overset{\rightarrow}{i}} + {\left( {\mathcal{G}_{4}^{1} + {\frac{\left( {\mathcal{G}_{2}^{1} - \mathcal{G}_{4}^{1}} \right)}{\Delta}x}} \right) \cdot \overset{\rightarrow}{j}}} \right\rbrack}}\quad\Delta\quad t}},} \\{{x \in \gamma_{p}},{t \in \left\lbrack {0,{\Delta\quad t}} \right\rbrack}}\end{matrix}$

In order to obtain the local function {tilde over (q)}(x, t) Equation(25) is applied:{tilde over (q)}(x, t)=λ(x, t){tilde over (g)}(x, t)where λ(x, t)=g(|∇I(x,t)|²). As ∇I(x,t) is a spatially sampled imagewhere samples are located at the 0-pixels of K^(p), the local values ofλ(x, t) are approximated. For the sake of simplicity, a bilinearapproximation is once again used; that is:{tilde over (λ)}(x)=a+bx+cy+dxy

For a 2-pixel of K^(p), as illustrated in FIG. 18, the followingrelation is obtained: $\begin{matrix}{{\overset{\sim}{\lambda}\left( {x,t} \right)} = {\lambda_{0,0} + {\frac{1}{\Delta}\left( {{\left( {\lambda_{1,0} - \lambda_{0,0}} \right)x} + {\left( {\lambda_{0,1} - \lambda_{0,0}} \right)y}} \right)} + {\frac{1}{\Delta^{2}}\left( {\lambda_{0,0} - \lambda_{1,0} - \lambda_{0,1} + \lambda_{1,1}} \right)}}} & (60)\end{matrix}$

Using these assumptions, the same steps as in hereinabove are followedto find Q as a function of G. For instance, this function for onen-pixel as defined in FIG. 16 c is:Q ₁=(CL)^(t) Gwhere C, L and G are matrices defined as: $\begin{matrix}{{C = \begin{bmatrix}{1/24} & {7/24} & {1/24} & {1/24} & {7/24} & {1/24} \\0 & {1/24} & {1/48} & 0 & {1/24} & {1/48} \\{1/48} & {1/24} & 0 & {1/48} & {1/24} & 0\end{bmatrix}},} \\\begin{matrix}{L = \begin{bmatrix}\lambda_{0,{- 1}} \\\lambda_{0,0} \\\lambda_{0,1} \\\lambda_{1,{- 1}} \\\lambda_{1,0} \\\lambda_{1,1}\end{bmatrix}} & \quad & {and} & \quad & {G = \begin{bmatrix}\mathcal{G}_{1}^{1} \\\mathcal{G}_{3}^{1} \\\mathcal{G}_{5}^{1}\end{bmatrix}}\end{matrix}\end{matrix}$

Using Equations (40) and (43), we can express D can be expressed as afunction of T. For one 3-pixel, γ_(E) of K^(s′) as defined in FIG. 19,and the following relation is obtained:<D, γ _(E)>=(C _(λ) L _(λ))^(t) TΔt   (61)where C, L and T are matrices defined as $\begin{matrix}{C_{\lambda} =} \\{\quad{\begin{bmatrix}{1/24} & {1/16} & 0 & {1/16} & {1/12} & 0 & 0 & 0 & 0 \\{1/48} & {1/2} & {1/48} & 0 & {5/24} & 0 & 0 & 0 & 0 \\0 & {1/16} & {1/24} & 0 & {1/12} & {1/16} & 0 & 0 & 0 \\{1/48} & 0 & 0 & {1/4} & {5/24} & 0 & {1/48} & 0 & 0 \\{{- 1}/12} & {{- 3}/8} & {{- 1}/12} & {{- 3}/8} & {{- 7}/6} & {{- 3}/8} & {{- 1}/12} & {{- 3}/8} & {{- 1}/12} \\0 & 0 & {1/48} & 0 & {5/24} & {1/4} & 0 & 0 & {1/48} \\0 & 0 & 0 & {1/16} & {1/12} & 0 & {1/24} & {1/16} & 0 \\0 & 0 & 0 & 0 & {5/24} & 0 & {1/48} & {1/4} & {1/48} \\0 & 0 & 0 & 0 & {1/12} & {1/16} & 0 & {1/16} & {1/24}\end{bmatrix},}} \\\begin{matrix}{L_{\lambda} = \begin{bmatrix}\lambda_{{- 1},{- 1}} \\\lambda_{0,{- 1}} \\\lambda_{1,{- 1}} \\\lambda_{{- 1},0} \\\lambda_{0,0} \\\lambda_{1,0} \\\lambda_{{- 1},1} \\\lambda_{0,1} \\\lambda_{1,1}\end{bmatrix}} & \quad & {and} & \quad & {T = \begin{bmatrix}\mathcal{T}_{{- 1},{- 1}}^{1} \\\mathcal{T}_{0,{- 1}}^{1} \\\mathcal{T}_{1,{- 1}}^{1} \\\mathcal{T}_{{- 1},0}^{1} \\\mathcal{T}_{0,0}^{1} \\\mathcal{T}_{1,0}^{1} \\\mathcal{T}_{{- 1},1}^{1} \\\mathcal{T}_{0,1}^{1} \\\mathcal{T}_{1,1}^{1}\end{bmatrix}}\end{matrix}\end{matrix}$

Equation 33 with <S, γ_(E)>=0 is applied to obtain, for each 3-pixel ofK^(p′):T _(0,0) ¹−(C _(λ) L _(λ))^(t) TΔt=T _(0,0) ⁰which defines the system of linear equations. The initial conditions areT⁰=I(x).A Different Hypothesis for Heat Conduction

In the preceding discussion, λ(x) has been approximated with a bilinearfunction, essentially for the sake of simplicity. Nevertheless, it couldbe preferable to use another assumption. Actually, this simple approachdoes not accurately handle abrupt changes in conductivity. For example,two 2-pixels of K^(s), as shown in FIG. 20 will be considered. Tocompute Q, λ(x) is approximated at the borders of the pixels based onthe values at their centers. Using bilinear approximation, the value of{tilde over (λ)}(x) on the line linking two points is declared be thearithmetic mean of the values at these points. For instance, givenλ_(0,0)→0 and λ_(1,0)→1, the conduction at the border is about 0.5. Thismeans that the zero conductivity at one pixel is partly cancelled out bythe fact that on the pixel beside it, there is a high conductivitycoefficient. In non-linear gray level diffusion, we are confronted withprecisely this kind of abrupt change. For example, at step edge pixels,the conduction may need to be very low, whereas immediately adjacent, itmay needs to be almost one.

A better assumption would thus be to consider {tilde over (λ)}(x) asconstant over one single 2-pixel of K^(s) [37]. Therefore, on the 1-facecommon to two pixels as in FIG. 20:${\overset{\sim}{\lambda}(x)} = {\left( {\frac{0.5}{\lambda_{0,0}} + \frac{0.5}{\lambda_{1,0}}} \right)^{- 1} = \frac{2\quad\lambda_{0,0}\lambda_{1,0}}{\lambda_{0,0} + \lambda_{1,0}}}$

It can easily be seen that when λ_(0,0)→0, then {tilde over (λ)}→0 andwhen λ_(0,0)<<λ_(1,0), then {tilde over (λ)}→2λ_(0,O). This means thatin both situations, the low conductivity would prevail at the boundarycommon to the two pixels [37]. With this assumption, the matrices C_(λ)and L_(λ) are modified as follows: $\begin{matrix}{{C_{\lambda} = \begin{bmatrix}{1/4} & {1/4} & 0 & 0 \\{3/2} & {{- 1}/4} & {{- 1}/4} & 0 \\{1/4} & 0 & {1/4} & 0 \\{{- 1}/4} & {3/2} & 0 & {{- 1}/4} \\{{- 3}/2} & {{- 3}/2} & {{- 3}/2} & {{- 3}/2} \\{{- 1}/4} & 0 & {3/2} & {{- 1}/4} \\0 & {1/4} & 0 & {1/4} \\0 & {{- 1}/4} & {{- 1}/4} & {3/2} \\0 & 0 & {1/4} & {1/4}\end{bmatrix}},{and}} \\{L_{\lambda} = {\lambda_{0,0}\begin{bmatrix}{\lambda_{0,{- 1}}/\left( {\lambda_{0,{- 1}} + \lambda_{0,0}} \right)} \\{\lambda_{{- 1},0}/\left( {\lambda_{{- 1},0} + \lambda_{0,0}} \right)} \\{\lambda_{1,0}/\left( {\lambda_{1,0} + \lambda_{0,0}} \right)} \\{\lambda_{0,1}/\left( {\lambda_{0,1} + \lambda_{0,0}} \right)}\end{bmatrix}}}\end{matrix}$Experimental Results

The proposed approach was tested on real and synthetic images in thecontext of linear isotropic diffusion, optical flow and non-lineardiffusion. The results were compared with another method in each case.

For linear diffusion, FIG. 21 a presents our physics-based method (PB)at three different scales and FIG. 21 b shows the result by convolutionfor the same scales. In the absence of a quantitative evaluation, it canbe said that subjectively the results seem to be similar.

For optical flow, FIGS. 22 a-22 c show the first frames of threesequences: rotating sphere, Hamburg taxi and tree sequences. The resultsare compared with those being obtained using a finite-differenceimplementation of the Horn and Schunck algorithm (FD) [18, 19]. In thesethree examples and for both the PB and FD methods, the image derivativesare computed by convolution with the appropriate Gaussian derivatives.Both temporal and spatial scales are set to 1, as is the weightingfactor α. FIGS. 23 a and 23 b shows the flow pattern computed for thesphere sequence. FIGS. 24 a-b and 25 a-b present the flow patterns forthe taxi and tree sequences, respectively. For the rotating sphere andthe taxi sequences, we obtain similar results with both methods. For thetree sequence, we also obtain similar results even if the extreme valuesseem to be smaller with the PB method than with the FD method. This factis more apparent in FIGS. 27 a and 27 b which show respectively theresults for the PB and FD methods for the tree sequence in which whitenoise (standard deviation of 10) has been added (see FIG. 26). Anotheradvantage of the method according to the present invention is that itavoids iterations since the algorithm is applied only once.

For nonlinear diffusion, FIGS. 28 a-28 c compare the PB method withconstant hypothesis on λ and the FD [38, 17] method for a small windowof the peppers image with σ=5. FIG. 28 b presents the original sectionwith white noise added (standard deviation of 10). FIGS. 28 b and 28 cshow respectively the results for the PB and FD methods. One can noticethat some details are better conserved in FIG. 28 c than in FIG. 28 b.This fact is enlightened in FIGS. 29 a-29 d which show a profile of thediagonal line starting at the upper right corner and finishing at thelower left corner.

FIGS. 30 a and 30 b present the results for the peppers image withσ=1.0, 3.0 and 5.0. The results for PB seem a little sharper than the FDresults.

FIGS. 32 a and 32 b show the results for the Lena image (FIG. 31 a) withan added white noise of standard deviation 10 (FIG. 31 b) at twodifferent scales, σ=4.0 and 8.0. Again, the PB method seems to givesharper results at both scales.

FIGS. 33 a and 33 b present details of the Lena results at σ=8.0. FIG.33 b seems smoother in constant zones but some details are lost. Forexample, compare the eyes, the trim on the hat and the right side of theface.

Conclusion

An alternative approach to the PDE-based resolution of the diffusionproblem was described. The proposed approach differs in two significantways from with the classical PDE resolution scheme: 1) the image isconsidered as a cubical complex for which algebraic structures such aschains, cochains, boundaries and coboundaries are defined; and 2) thediffusion problem is decomposed into conservative and constitutive basiclaws, each of which is represented by cochains and coboundaries.

The conservative basic laws are represented without approximation whilesome approximations are required for the constitutive laws. This meansthat unlike traditional PDE resolution, for which many approximationsmust be made, all approximations are known since they are only needed inthe representation of the constitutive equations. Coboundaries arecomputed using fundamental theorems of calculus such as the Green,Stokes and Line Integral theorems. Unlike iterative numerical analysisalgorithms that do not allow the explanation of intermediate results,the use of basic laws allows the physical explanation of all steps andintermediate results of the algorithm. Moreover, since there is noiteration in the resolution process, there is no problem about theconvergence of the numerical analysis scheme. Furthermore, the use ofcubical complexes provides algorithms that can operate in any dimension.It has the significant advantage of avoiding the potentially difficulttask of extending the algorithm to higher dimensions. Cochains andcoboundaries allow the use of both global and local quantities.Integrals or discrete summations over fields are used to compute globalquantities. This allows the reduction of noise by performing a smoothingoperation, as opposed to differentiation, which enhances highfrequencies.

In computer vision and image processing, several problems can be modeledas diffusion problems. The proposed approach has been validated onsmoothing by linear and nonlinear diffusion and on the computation ofoptical flow. The results obtained confirm the effectiveness of thisapproach.

PRACTICAL EXAMPLE #2 A Physics-Based Model for Active Contours

A new active contours model is presented. It is based upon adecomposition or the linear elasticity problem into basic physical laws.As opposite with other physics-based active contours model which solvethe partial differential equation arising from the physical laws by somepurely numerical techniques, exact global values are used andapproximations made only when they are needed. Moreover, theseapproximations can be made wisely assuming some knowledge about theproblem and the domain. The deformations computed with the presentapproach have a physical interpretation. In addition, the deformedcurves have some interesting physical properties such as the ability torecover their original shape when the external forces are removed. Thephysical laws are encoded using the computational algebraic topologybased image model described herein. The resultant numerical scheme isthen straightforward. The image model allows our algorithm to performwith either 2D or 3D problems.

Introduction

These last years, active contours and active surfaces have been widelystudied since the introduction of active contours by Kass et al [59].They have been used in image segmentation [62], tracking [68], automaticcorrection and updating of road databases [46], etc.

To solve these problems, many different approaches have been proposed(see [57, 63]) in particular physical models derived from equations ofcontinuum mechanics. Mass-springs models are physical models which use adiscrete representation of the objects. Objects are modeled as a latticewith point masses linked together by springs [57]. Information is thusonly available at a finite number of points [63]. These methods offeronly a rough approximation of the phenomena happening in a body [57].Moreover, the determination of spring constants reflecting the materialproperties may be a very fastidious work. However, they offer real-timeperformances and allows for parallel computations.

Other physical models based upon the minimization of an energyfunctional which takes into account an internal regularizing force andan external force applied on the data are also often used. Some of themconsider the deformable bodies as continuous objects by approximatingtheir continuous behavior with methods such as the finite element method(FEM). FEM are closer to the physics than mass-springs models but theircomputational requirements make them difficult to be applied inreal-time systems without preprocessing steps [57]. Finite differencemethods (FDM) are also used to discretize the objects. They usuallyoffer better performance than FEM but they require the computation offourth order derivatives which make them sensitive to noise [63].

For a given curve S₁ the application of the FEM and FDM methods leads toa discrete stationary system of equations of the form:KS=f(S)where K is a matrix which encodes the regularizing constraints on S andf(S) represents the data potential. However, some problems such asanimation in graphics applications require to take into account adynamic evolution of the curve [57]. In these case, inertial body forcesand damping forces may also by considered by controlling thedeformations through a Newtonian law of motion: $\begin{matrix}{{{M\quad\frac{\partial^{2}S}{\partial t^{2}}} + {D\quad\frac{\partial S}{\partial t}} + {K\quad S}} = {f(S)}} & (62)\end{matrix}$or by a Lagrangian evolution $\begin{matrix}{{\frac{\partial S}{\partial t} + {K\quad S}} = {f(S)}} & (63)\end{matrix}$where M and D are respectively matrices which represent the mass modeland the background damping. Equations (62) and (63) are solved usingvarious numerical schemes [70, 64] assuming an initial curve S₀ close tothe solution which evolves until the inertial terms go to zero.

Over the last years, a lot of different methods have been introduced tocompute the matrices M, D and K but as pointed out by Montagnat et al[63], these methods have a major drawback: the corresponding systemdeformations do not have any physical interpretations.

A new model which includes a physical interpretation of the deformationsis described hereinbelow. The model is similar to a mass-springs modelbut it includes both the efficiency of the mass-springs models and theaccuracy of the physical modeling of the FEM by providing a systematicmethod for specifying springs constants reflecting the properties of thematerials.

To achieve it, we propose to use directly the basic laws of physicswhich lead to the partial differential equations (62) and (63). Theseequations are indeed obtained by considering and mixing together somebasic laws of physics into a global conservation law and considering itslocal counterpart [72]. This approach is not always well suited forproblems such as computer vision in which the continuous domain must besubdivided into many sub-domains for which there are often only oneinformation available. The use of this information as a global valueover each sub-domain allow to directly use the global conservation lawwhich can lead to an algorithm less sensitive to noise.

To encode these global values over points, surfaces, volumes, etcarising from some physical laws, we use the computational algebraictopology based image model described hereinabove.

The approach according to this illustrative embodiment of the presentinvention has several advantages. 1) Since the linear elasticity problemis well-known in continuum mechanics, the modeling according to theillustrative embodiment of the present invention can be made wisely inorder to provide some good physical interpretation of the wholedeformation process and of its intermediate steps. This allows an easierdetermination of the parameters used in the process since they have aphysical meaning; 2) The determination of the springs constants in orderto reflect the material properties is straightforward; 3) The objects inthe image (e.g. curves, surfaces) are modeled as entities having theirown physical properties such as elasticity and rigidity. They have theproperty of recovering their original state when the forces applied onthem are removed; 4) Both smooth results and results having highcurvature points can be obtained; 5) The complexity of the algorithm isminimal and allows for real-time simulation without any extrapreprocessing steps [51]; 6) The image model allows our algorithm toperform with either 2D or 3D problems.

Physical Modeling

One of the objectives of this illustrative embodiment of the presentinvention is to model the objects in an image (e.g. curves, surfaces,etc) as entities having their own physical properties such as elasticityand rigidity. As a consequence, these objects need to satisfy the lawsand principles of the continuum mechanics. For instance, a bodysubjected to forces must move or deform according to the universal lawsof physics.

These principles and laws to which all bodies must obey will now bepresented. We first introduce the concepts of stress and strain whichare required in the statement of the governing equations for deformablebodies. Then, the physical laws related to the linear elasticity problemwill be presented.

The elasticity theory has been widely studied by engineers andscientists and is the main subject of many books. The presentspecification only presents the concepts of this theory which arerelevant to our application. The concepts presented here are well knownand may be found in many continuum mechanics books such as [65, 50].

Forces, Stresses and Strains

A material body in a 3-D space is always subjected to forces. Theseforces may come from an external agent (external forces) or issue fromthe object itself (internal forces). When the external forces aregreater than the internal forces, then the body can undergo deformations(strains) or be accelerated. This deformation can induce internal forces(stresses) if the material is elastic. The concepts of force, stress andstrain and the relation between strain and stress is exposedhereinbelow.

Forces Acting on a Body

Two basic types of forces act on a body. First, there are theinteratomic forces which hold the body's particles together at someconfiguration. These forces, called internal forces, can either attemptto separate or bring the particles closer according to the fact that thebody undergoes a contraction or an extension [65]. They act in responseto a force applied by some other agent. Assuming the equilibrium of thebody and using Newton's law of reaction, they must be equal in magnitudeto the forces applied to the body but in opposite direction [66].

On the other hand, there are the forces applied by an external agent,called external forces. Two types of these forces are generally appliedon a body. First, there are forces such as gravity and inertia calledthe body forces, which act on all volume elements. These forces, notedb_(i) (forces per unit of mass in a direction x_(i)), are distributed inevery part of the body. Secondly, there are forces which act on and aredistributed over a surface element such as the contact forces betweensolid elements [49]. These forces, noted fi (forces per unit of area ina direction x_(i), are called the surface forces. The surface elementmay be inside the body or a part of a bounding surface [60]. A body ofarbitrary size, shape and material subjected to surface and body forcesis shown in FIG. 34.

External forces applied on a body must be transmitted to it. A rigidbody can then undergo either a spatial shift, a rotation of both ofthem. In the case of a non-rigid body, it can also go through adeformation or a distortion in which case internal forces will bedeveloped to counterbalance the external forces. If the internal andexternal forces are balanced, we say that the body is in staticequilibrium. Otherwise the body can be accelerated which would give riseto inertia forces [58]. Using d'Alembert's principle, these forces maybe included as part of the body forces [48] such that the equilibriumequations can be satisfied. If the body gets deformed then thedeformation can either be elastic or not and is subjected to thematerial properties of the body such as its elasticity and its rigidity.If the internal forces induced by the material properties of a body aretoo weak to counterbalance the external forces, then the body can bepermanently deformed [52].

We assume herein the material to be isotropic with respect to somemechanical properties. We then suppose that the material properties arethe same in all directions for a given point [60]. We also consider anhomogeneous material which means that its properties are identical atall locations.

The Concept of Stress at a Point

Let us consider an isotropic and homogeneous body B. Let us assume thatB is subjected to arbitrary surface and body forces such that B is instatic equilibrium. Let P be a interior point of B and S be a planesurface passing through P. S will be referred to as the cutting planeand is defined by the unit normal vector n=(n₁, n₂, n₃)^(T). Then Spartitions B into two sections I and II as shown in FIG. 35.

Let us assume that ΔS is a small element of area of the cutting planesurrounding P (see FIG. 36).

Since the body is in static equilibrium then the force system acting oneach part I and II taken alone must also be in equilibrium. Thisgenerally requires that some internal forces are transmitted by part Ito part II. These forces are not necessarily distributed uniformly onevery part of the cutting plane. Thus they may vary in magnitude anddirection over it. We generally want to determine precisely that forcedistribution at every point of ΔS. The term stress is used to define theintensity and the direction of the internal force Δf acting at point P.Using the Cauchy stress principle [60] we define the stress vector (ortraction vector or traction forces) t^(n) at P as:$t^{n} = {\lim\limits_{{\Delta\quad S}->0}\frac{\Delta\quad f}{\Delta\quad S}}$assuming that P remains an interior point of ΔS as its area reduces tozero.

Let us mention that t^(n) is not necessarily in the direction of thenormal vector n at P. However, it may be decomposed into a componentperpendicular to the cutting plane, called the normal stress, and acomponent parallel to it, called the shear stress. The normal stressattempts to separate (bring closer) the material particles after acompression (an extension) of an elastic body when it tries to recoverits original state. On the other hand, the shear stress acts parallel tothe cutting plane and tends to slide adjacent planes with respect toeach other (see FIG. 37 a-37 d).

It should be first noticed that the stress vector t^(n) is defined withrespect to the cutting plane's unit normal vector n. Since there areinfinitely many cutting planes going through P there are also as manystress vectors defined at P. Juvinall [58] defines the state of stressat a point P as a complete description of the stress magnitude anddirection for all possible cutting plane passing through P. Fortunately,this description can be fully obtained by considering any three mutuallyperpendicular planes passing at P [58]. For the sake of simplicity, weusually use the three axes defined by the three canonical vectors x₁, x₂and x₃.

Let us define σ_(ij) as the stress component in the direction of x_(j)when the normal vector is parallel to the axe defined by x_(i). If i=jthen σ_(ii) represents a normal stress. Otherwise, σ_(ij) is a shearstress. With these conventions, the component t_(i) ^(n) in thedirection of x_(i) of the traction force t^(n) depends on the normalstress σ_(ii), the shear stresses σ_(ji) and σ_(ki) and the normalvector n such that (see FIG. 38): $\begin{matrix}\begin{matrix}{t_{i}^{n} = {{\sigma_{1i}n_{1}} + {\sigma_{2i}n_{2}} + {\sigma_{3i}n_{3}}}} \\{= {\sum\limits_{j = 1}^{3}{\sigma_{ji}n_{j}}}}\end{matrix} & (64)\end{matrix}$Equation (64) is known as the Cauchy stress formula.

Since each of the three coordinates axes involves six stress components,there is a total of nine stress components. However the equilibrium ofmoments at P [49] gives that only six of these are independent that isσ_(ij)=σ_(ji) for all i, j=1, 2, 3. This means that the state of stressat a point is fully determined by σ₁₁, σ₂₂, σ₃₃, σ₁₂, σ₁₃ and σ₂₃.

The Concept of Strain at a Point

Any non-rigid body goes through deformations and distortions whensubjected to forces. The body can either extend or contract(deformation) or have a geometric modification of its shape(distortion). FIGS. 39 a and 39 b present these concepts.

The term strain refers to the direction and intensity of the deformationat any given point with respect to a specific plane passing through thatpoint [58]. As for stress, the strain is defined according to a specificcutting plane. The state of strain is defined by Juvinall [58] as thecomplete definition of the magnitude and direction of the deformation ata given point with respect to a all cutting planes passing through thatpoint.

As for the state of stress, the description can be obtained byconsidering any three mutually perpendicular planes passing at P. Onecan therefore see a great similarity between stress and strain. However,there is a major difference between them: strains are generally somedirectly measurable quantities while stresses are not. Fortunately,stresses can be computed from strains (and vice versa) using aconstitutive equation.

As for stress, two types of strains can be defined. First, there are thestrains which result of a change in the dimensions of the body(deformation). Let B be the same body defined hereinabove, ΔB be a smallelement of B of length Δx_(i) in the direction of x_(i) and ΔB′ be thedeformation of ΔB such that Δu_(i) is the change in lenght of ΔB afterthe application of a force in the direction of x_(i) (see FIG. 40). Thenormal strain ε_(il) at P in the x_(i) direction with respect to acutting plane having x_(i) as normal vector is the unit deformation of aline element [65] in the direction of x_(i). It is formally defined as:$ɛ_{ii} = {{\lim\limits_{{\Delta\quad x_{i}}->0}\frac{\Delta\quad u_{i}}{\Delta\quad x_{i}}} = \frac{\partial u_{i}}{\partial x_{i}}}$

The normal strain is clearly the unit change in length per unit originallength for the element in the direction of x_(i). Since it is the ratioof two units of length, it is dimensionless even if it is sometimesexpressed as units of length per unit of length such as inches per inch.

Let us now suppose that there are two perpendicular lines PB and PA oflength Δx_(j) and Δx_(k) respectively in the direction of x_(j) andx_(k) (see FIG. 41).

Let us assume that after a distortion points A and B move respectivelyto A′ and B′ while P remains fixed. The lines PA and PB have beenrotated of angles θ_(jk) and θ_(kj) such that:$\frac{\Delta\quad u_{j}}{\Delta\quad x_{k}} = {{{\tan\left( \theta_{jk} \right)}\quad{and}\quad\frac{\Delta\quad u_{k}}{\Delta\quad x_{j}}} = {\tan\left( \theta_{kj} \right)}}$where Δu_(j) and Δu_(k) are respectively the displacements of B and A inthe x_(j) and x_(k) directions.

If it is assumed that only small distortions occur, then we canapproximate both tangents by their angles. Thus: $\begin{matrix}{\theta_{jk} \cong {\tan\left( \theta_{jk} \right)}} \\{= {\frac{\Delta\quad u_{j}}{\Delta\quad x_{k}}\quad{and}}} \\{\theta_{kj} \cong {\tan\left( \theta_{kj} \right)}} \\{= \frac{\Delta\quad u_{k}}{\Delta\quad x_{j}}}\end{matrix}$or, taking their infinitesimal analogous: $\begin{matrix}\begin{matrix}{\theta_{jk} = {\lim\limits_{{\Delta\quad x_{k}}->0}\frac{\Delta\quad u_{j}}{\Delta\quad x_{k}}}} \\{= {\frac{\partial u_{j}}{\partial x_{k}}\quad{and}}} \\{\theta_{kj} = {\lim\limits_{{\Delta\quad x_{j}}->0}\frac{\Delta\quad u_{k}}{\Delta\quad x_{j}}}} \\{= \frac{\partial u_{k}}{\partial x_{j}}}\end{matrix} & (65)\end{matrix}$

The shear strain γ_(ik) at P with respect to the cutting plane havingx_(i) as normal vector is the angle in radians through which twoorthogonal lines in the undistorted body are rotated by a distortion[56]. That is γ_(ik)=θ_(jk)+θ_(kj). The two subscripts in γ_(ik) have asimilar meaning as for stress. For instance, γ_(ik) is the strain actingon two adjacent planes perpendicular to the x_(i) axis and sliding themrelative to each other in the x_(k) direction.

Let us recall that these definitions have been made under the assumptionthat only very small displacements occur in the body. The normal andshear strains are supposed small compared to unity [50]. If thisconstraint is relaxed in order to include large deformations then thesystem to solve for the computation of the forces, the stresses, thestrains or the displacements becomes non-linear and then harder tosolve. This is sometimes necessary in some problems where largedeformations can occur such as for thin flexible bodies [50] of for themodelization of human tissue [53]. However, if we restrict ourselves tosmall deformations, then the approximations made to define the shearstrains do not induce too many errors and are widely accepted in theclassical theory of elasticity [69, 65, 49, 50].

Finally, Equation (65) clearly shows that the shear strain is also adimensionless quantity since it is the ratio of two units of length.Defining for i≠j:$ɛ_{ij} = {\frac{1}{2}\gamma_{ij}\quad\left( {i,{j = 1},2,3} \right)}$the followin strain-displacement relation (or the kinematicalrelationship [69]) is obtained: $\begin{matrix}{ɛ_{ij} = {\frac{1}{2}\left\lbrack {\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}}} \right\rbrack}} & (66)\end{matrix}$

As for stresses, there is a total of nine strain components at eachpoint of the body (three per mutually perpendicular cutting planes) butby symmetry they can be reduced to six, that is γ_(jk)=γ_(kj) for all j,k=1, 2, 3 with j≠k. The state of strain at any point can then bedescribed by ε₁₁, ε₂₂, ε₃₃, ε₁₂, ε₁₃ and ε₂₃.

Relations Between Forces, Stresses and Strains

As mentioned hereinabove, strains are measurable quantities while stressare not even if both can be computed from the other. This is due to thefact that some knowledge about the material of a body is necessary tomeasure the stresses from strains and vice versa. For instance, a steelbeam and a rubber beam induce different internal forces when bentequally.

This gap between strains and stresses will now be filled by stating thephysical laws relating them to each other. This hole between strains andstresses needs to be filled by a constitutive equation (or materiallaw), which reflect the internal constitution of the materials. Thematerial law for the linear elasticity problem is known as the Hooke'slaw.

Before stating the law, it should be first reminded that a material hasan elastic behavior when it satisfies the two following conditions:

-   -   1. The stresses depend only on the strains.    -   2. Its properties allow a body to recover its original shape        when the external forces applied on the body are removed [60].

If theses conditions are not satisfied, a body is said to have aninelastic behavior. Any body may be seen as having an elastic behavioras long as it is not deformed beyond a limit value. This value is calledthe elastic limit [52, 49] and is usually defined as the maximum valueof stress that a body can undergo without undergoing a permanentdeformation.

In addition, if the stress is a linear function of the strain, anelastic material has a linear elastic behavior. In what follows, it isassumed that the material has a linear elastic behavior.

As mentioned in [60] and [49], in many situations the problem ofelasticity can be considered as a 2D problem. This particular case isknown as plane elasticity. Two basic types of problems compose the planeelasticity. The problems in which the stress components in one directionfor a body are all zero are referred to as plane stress problems. On theother hand, if all the strain components in one direction for a body arezero then the state of strain for that body is referred to as planestrain (see [65, 60, 58]).

The present problem is considered as a plane strain problem. Thisdistinction is important since the constitutive equation slightlydiffers according to the fact that a plane stress or a plane strainproblem is considered.

The Hooke's Law

When a rubber ball is compressed its diameter in the directionsperpendicular to the applied force gets larger. A similar phenomenonoccurs when a rubber band is extended and its cross section getssmaller¹. In fact, these changes in dimensions happen in all materialseven if they can't always be noticed by a naked eye [52].¹ Example taken in [52]

When a stress is acting on an isotropic and homogeneous body in only onedirection (uniaxial stress), one can show that the transverse strainε^(⊥) (the strain in a perpendicular direction) is directly proportionalto the strain ε induced by the stress:ε^(⊥) =−vεThe ratio: $v = \frac{- ɛ^{\bot}}{ɛ}$is called the Poisson ratio. It is supposed constant when the stress isbelow the elastic limit [66].

The linear relationship between a uniaxial stress σ_(ii) in a directionx_(i) and the corresponding strain ε_(ii) is known as the Hooke's law[52, 58] and is written as: $ɛ_{ii} = \frac{\sigma_{ii}}{E}$where E is the Young's modulus of elasticity. Of course, the Hooke's lawis valid if the stress is not beyond the elastic limit of the material.

As pointed out by Boresi [49], the stresses at any point depend on allthe strains in the neighborhood of that point. Thus the totaldeformation in the x_(i) direction depends not only on the stress inthat direction but also of the deformations in the other twoperpendicular directions. For instance the normal strain ε_(ii) does notonly depend on $\frac{\sigma_{ii}}{E}$(Hooke's law) but also of the transverse strains ε_(jj) and ε_(kk) suchthat the total deformation in the direction of x_(i) is: $\begin{matrix}\begin{matrix}{ɛ_{ii} = {\frac{\sigma_{ii}}{E} - {v\quad ɛ_{jj}} - {v\quad ɛ_{kk}}}} \\{= {\frac{\sigma_{ii}}{E} - {v\quad\frac{\sigma_{jj}}{E}} - {v\quad\frac{\sigma_{kk}}{E}}}} \\{= {{\frac{1}{E}\left\lbrack {\sigma_{ii} - {v\left( {\sigma_{jj} + \sigma_{kk}} \right)}} \right\rbrack}\quad\left( {i \neq j \neq k} \right)}}\end{matrix} & (67)\end{matrix}$

Equation (67) is the normal strain-stress relation. For isotropicmaterials, it can be shown that the normal strains are not influenced bythe shear stresses [52]. Consequently, the shear stresses only induceshear strains and they are related by the relation: $\begin{matrix}{{2ɛ_{ij}} = {\frac{\sigma_{ij}}{G}\quad\left( {i \neq j} \right)}} & (68)\end{matrix}$where $G = \frac{E}{2\left( {1 + v} \right)}$is called the modulus of rigidity. Equations (67) and (68) are known asthe Generalized Hooke's law for linear elastic isotropic materials [52].These equations can be inverted in order to express the stresses asfunctions of the strains: $\begin{matrix}{\sigma_{ii} = {\frac{E}{\left( {1 + v} \right)\left( {1 - {2v}} \right)}\left\lbrack {{\left( {1 - v} \right)ɛ_{ii}} + {v\left( {ɛ_{jj} + ɛ_{kk}} \right)}} \right\rbrack}} & (69) \\{\sigma_{ij} = {{2G\quad ɛ_{ij}} = {\frac{E}{1 + v}ɛ_{ij}\quad\left( {i \neq j} \right)}}} & (70)\end{matrix}$Relation Between Forces and Stresses

It is well known that the conservation (balance, equilibrium) lawsconstitute an important class of equations in continuum mechanics. Theyrelate the change in total amount of a physical quantity inside a bodywith the amount of this quantity, which flows through its boundary.These laws must be satisfied for every continuous materials. Localdifferential equations are usually used to express these laws. In whatfollows, the linear momentum, which is relevant for the linearelasticity problem, is presented.

To every material body B is associated a measure of its inertia calledthe mass. This measure may vary in space and time inside a body. Let Vbe the volume of B, S its bounding surface and Δm be the mass of a smallamount of volume ΔV. The mass density is given by:$\rho = {{\rho\left( {x,t} \right)} = {\lim\limits_{{\Delta\quad V}->0}\frac{\Delta\quad m}{\Delta\quad V}}}$

Let us assume that distributed body forces ρb_(i) and tractions forcest_(i) ^(n) are applied to S (see FIG. 42). Let also assume that B ismoving under the velocity field v_(i)=v_(i)(x,t). The quantity:${P_{i}(t)} = {\underset{V}{\int{\int\int}}\rho\quad\upsilon_{i}{\mathbb{d}V}}$is called the linear momentum of B. The principle of linear momentum[49] states that the resultant force acting on a body is equal to thetime rate of change of the linear momentum. Thus: $\begin{matrix}{{\frac{\mathbb{d}}{\mathbb{d}t}\underset{V}{\int{\int\int}}\rho\quad\upsilon_{i}{\mathbb{d}V}} = \underset{\underset{{Forces}\quad{acting}\quad{on}\quad{the}\quad{body}}{︸}}{{\underset{S}{\int\int}t_{i}^{n}{\mathbb{d}S}} + {\underset{V}{\int{\int\int}}\rho\quad b_{i}{\mathbb{d}V}}}} & (71)\end{matrix}$Recalling Equation (64):$t_{i}^{n} = {\sum\limits_{j = 1}^{3}{\sigma_{ji}n_{j}}}$where n=(n₁, n₂, n₃)^(T) is the unit normal vector to the surface andσ=(σ₁₁, σ₂₁, σ₃₁)^(T) and using Gauss's divergence theorem, thefollowing relation is obtained: $\begin{matrix}{{\frac{\mathbb{d}}{\mathbb{d}t}\underset{V}{\int{\int\int}}\rho\quad\upsilon_{i}{\mathbb{d}V}} = {{\underset{V}{\int{\int\int}}\rho\quad\frac{\partial\upsilon_{i}}{\partial t}{\mathbb{d}V}} = {\underset{V}{\int{\int\int}}\left( {{\nabla{\cdot \sigma}} + {\rho\quad b_{i}}} \right){\mathbb{d}V}}}} & (72)\end{matrix}$

Since V is arbitrary then the integral sign can be retrieved leading tothe local equations of motion: $\begin{matrix}{{\rho\frac{\mathbb{d}\upsilon_{i}}{\mathbb{d}t}} = {{\nabla{\cdot \sigma}} + {\rho\quad b_{i}}}} & (73)\end{matrix}$

The global equilibrium equations can be obtained assuming a zerovelocity field in Equation (72): $\begin{matrix}{{\underset{\underset{{Internal}\quad{forces}}{︸}}{\int{\int_{V}{\int{{\nabla{\cdot \sigma}}{\mathbb{d}V}}}}} + \underset{\underset{{External}\quad{forces}}{︸}}{\int{\int_{V}{\int{\rho\quad b_{i}{\mathbb{d}V}}}}}} = {0\quad\left( {{i = 1},2,3} \right)}} & (74)\end{matrix}$and their local counterparts are:∇·σ+ρb _(i)=0 (i=1, 2, 3)   (75)

Let us now summarize the decomposition of the problem. The relationbetween the global displacement U of a body and the correspondingstrains (see FIG. 43 a) has been introduced hereinabove. Theconstitutive equation relating the strains and the stresses (see FIG. 43b) has also been presented. Finally, how the stresses are related toforces using the linear momentum principle (see FIG. 43 c) has beendescribed. A general scheme similar to the one presented by Tonti [72]may then be introduced to summarize how the internal reaction forces ofa body are related to the global displacements of that body (FIG. 44).

Discrete Representation of Images

Some algebraic tools used to model the image will now be recalled fromthe above description. An image is composed of two distinctive parts:the image support (pixels) and some field quantity associated to eachpixel. This quantity can be scalar (e.g. gray level), vectorial (e.g.color, multispectral, optical flow) or tensorial (e.g. Hessian). Theimage support is modelled in terms of cubical complexes, chains andboundaries. With these concepts, it is possible to give a formaldescription of an image support of any dimension. For quantities, theconcept of cochains which are representations of fields over a cubicalcomplex is introduced. For the use of these concepts in imageprocessing, see [45].

An image is a complex of unit cubes usually called pixels. A pixelγ⊂R^(n) n is a product:γ=I ₁ ×I ₂ × . . . ×I _(n)where I_(j) is either a singleton or an interval of unit length withinteger end points. Then I_(j) is either the singleton {k} and is saidto be a degenerate interval, or the closed interval [k, k+1] for some kε Z. The number q ε {0,1, . . . , n} of non-degenerate intervals is bydefinition the dimension of γ which is called a q-pixel. FIGS. 45 a-45 billustrate three elementary pixels in R².

For q≧1, let J={k₀, k₁, . . . , k_(q−1),} be the ordered subset {1, 2, .. . , n} of indices for which I_(k) _(j) =[a_(j), b_(j)] isnon-degenerate. Define:A _(k) _(j) σ=I ₁ × . . . I _(k) _(j) ⁻¹ ×{a _(j) }×I _(k) _(j) ₊₁ × . .. ×I _(n)andB _(k) _(j) σ=I ₁ × . . . I _(k) _(j) ⁻¹ ×{b _(j) }×I _(k) _(j) ₊₁ × . .. ×I _(n)

The A_(k) _(j) and the B_(k) _(j) are called the (q−1)-faces of σ. Onecan define the (q−2)-faces, . . . , down to the 0-faces of σ the sameway. The faces of γ different from γ itself are called its proper faces.

By definition, a natural orientation of the cube is assumed for eachpixel. Suppose that γ denotes a particular positive oriented q-pixel. Itis natural to denote the same pixel with opposite orientation by −γ.Examples of orientations are given in FIGS. 45 a-45 c. A cubical complexin R^(n) is a finite collection K of q-pixels such that every face ofany pixel of the image support called K is also a pixel in K and theintersection of any two pixels of K is either empty or a face of each ofthem. For example, traditional 2D image models was only consideringpixel as a 2D squared element. Definitions presented before allows us toconsider 2-pixels (squared elements), 1-pixels (line elements) and0-pixels (punctual elements) simultaneously.

In order to write the image support in algebraic form, the concept ofchains is introduced. Any set of oriented q-pixels of a cubical complexcan be written in algebraic form by attributing them the coefficient 0,1or −1, if they are not in the set or if they should or not be taken withpositive orientation, respectively. In order to represent weighteddomains, arbitrary integer multiplicity for each q-pixel is allowed.

Given a topological space X ⊂ R^(n) in terms of a cubical complex, afree abelian group C_(q)(X) generated by all the q-pixels is obtained.The elements of this group are called q-chains and they are formallinear combinations of q-pixels [45]. A formal expression for a q-chainc_(q) is$c_{q} = {\sum\limits_{\gamma_{i} \in K}{\lambda_{i}\gamma_{i}}}$where λ_(i) ε Z.

The last step needed for the description of the image plan is theintroduction of the concept of boundary of a chain. Given a q-pixel γ,its boundary ∂γ is defined as the (q−1)-chain corresponding to thealternating sum of its (q−1)faces. The sum is taken according to theorientation of the (q−1)-faces with respect to the orientation of theq-pixel. It is said that a (q−1)-face of γ is relatively positivelyoriented with respect to the orientation of γ if its orientation iscompatible with the orientation of γ. By linearity, the extension of thedefinition of boundary to arbitrary q-chains is expedient. For instance,in

FIGS. 45 b and 45 c, the boundary of the 1-pixel a is x₂−x₁ and theboundary of the 2-pixel A is a+b−c−d whereby a and b are positivelyoriented with respect to orientation of A but c and d are negativelyoriented with respect to orientation of A. Let us notice that theboundary of a 1-pixel is always the difference of its boundary points.The boundary can be defined recursively. Suppose a (q−1)-chain and aq-chain γ_(q) defined as γ_(q)=γ_(q−1)×[a, b], the boundary of γ_(q) canbe recursively written as:∂γ_(q)=∂γ_(q−1) ×[a, b]+(−1)^((q−1))(γ_(q−1) ×{b{−γ_(q−1) ×{a})   (76)

In order to model the pixels quantity over the image plane, anapplication F has to be found to associate a global quantity with allq-pixels γ of a cubical complex. This application is denoted <F, γ>.This quantity may be any mathematical entities such as scalars, vectors,etc. For two adjacent q-pixels γ₁ and γ₂. F must satisfy <F,λ₁γ₁+λ₂γ₂>=λ₁<F, γ₁>+λ₂<F, γ₂>, which means that the sum of the quantityover each pixel is equal to the quantity over the two pixels. Theresulting transformation F:C_(q)(X)→R is called a q-cochain and is usedas a representation of a quantity over the cubical complex.

An operator is finally needed to associate a global quantity to the(q+1)-pixels according to the global quantities given on their q-faces.Given a q-cochain F, an operator ∂, called the coboundary operator, isused to transform F into a (q+1)-cochain ∂F such that:<δF, γ>=<F, ∂γ>  (77)for all (q+1)-chains γ. The coboundary is defined as the signed sum ofthe physical quantities associated with the q-faces of γ. The sum istaken according to the relative orientation of the q-faces of the(q+1)-pixels of γ with respect to their orientation. FIG. 46 presents anexample of the coboundary operation for a 2-pixel.Representation of the Equilibrium Equation

The basic laws of FIG. 44 with concepts of algebraic topology in orderto get a generic algorithm for solving the equilibrium Equation (74)will now be modeled.

The algorithm is resumed as follows: 1) The image support is firstlysubdivided into cubical complexes; 2) Global quantities are computedover pixels of various dimensions via cochains according to basic laws;3) The constitutive equations 69 and 70 are expressed as a lineartransformation between two cochains.

The Relative Displacement

Let B be a body in a 3D space and K^(p) be a 3-complex representing thesubdivided spatial support of B. Let us consider a 0-cochain U and a1-cochain D such that D is the coboundary of U: $\begin{matrix}\begin{matrix}{{\mathcal{D}:{\mathcal{C}_{1}\left( \mathcal{K}^{p} \right)}}->{\mathbb{R}}^{3}} \\{\left. \gamma\mapsto\left\langle {\mathcal{D},\gamma} \right\rangle \right. = {\left\langle {{\delta\mathcal{U}},\gamma} \right\rangle = \left\langle {\mathcal{U},{\partial\gamma}} \right\rangle}}\end{matrix} & (78)\end{matrix}$

FIGS. 47 a and 47 b present some examples of U and D for a 3-pixel ofK^(p).

The computational rules for both cochains U and D will now be specified.Recalling the strain-displacement relation (Equation (66)):$\begin{matrix}{ɛ_{ij} = {\frac{1}{2}\left\lbrack {\frac{\partial u_{j}}{\partial x_{i}} + \frac{\partial u_{i}}{\partial x_{j}}} \right\rbrack}} & (79)\end{matrix}$we have an application ε′: $\begin{matrix}{{ɛ^{\prime}:{\mathbb{R}}^{3}}->{\mathbb{R}}^{6}} \\{\left. U\mapsto{ɛ^{\prime}(U)} \right. = \left( {ɛ_{11},ɛ_{22},ɛ_{33},ɛ_{12},ɛ_{13},ɛ_{23}} \right)^{T}}\end{matrix}$

Omitting the shear strain components as in [67], the following relationmay be defined: $\begin{matrix}\begin{matrix}{{ɛ^{\prime}:{\mathbb{R}}^{3}}->{\mathbb{R}}^{3}} \\{\left. U\mapsto{ɛ(U)} \right. = {\left( {ɛ_{11},ɛ_{22},ɛ_{33}} \right)^{T} = {\nabla U}}}\end{matrix} & (80)\end{matrix}$

Using the global form of Equation (80) over a 1-pixel γ_(D) such that∂γ_(D)=x_(*)−x_(#), the following relation is obtained: $\begin{matrix}{{\int_{\gamma\quad D}{{ɛ(U)}{\mathbb{d}\gamma_{D}}}} = {{\int_{x_{\#}}^{x_{*}}{{ɛ(U)}{\mathbb{d}\gamma_{D}}}} = {\int_{x_{\#}}^{x_{*}}{{\nabla U}{\mathbb{d}\gamma_{D}}}}}} & (81)\end{matrix}$where dγ_(D) is an infinitesimal part of the domain γ_(D). Since Δu is aconservative field, applied is the line integral theorem [55, 71] whichstates that for a conservative field F(x)=∇f(x) and for two points A andB in an opened connected region containing F(x) the integral of thetangential part of F(x) along a curve R joining A and B Is independentof the path: ∫_(A)^(B)F(x)𝕕R = f(B) − f(A)

From Equation (81), the following relation is then obtained:$\begin{matrix}{{\int_{x_{\#}}^{x_{*}}{{\nabla U}{\mathbb{d}\gamma_{D}}}} = {{U\left( x_{*} \right)} - {U\left( x_{\#} \right)}}} & (82)\end{matrix}$

On the other hand, applying the cochain D to the 1-pixel γ_(D) leads to:<D, γ _(D) >=<U, ∂γ _(D) >=U(x _(*) −x _(#))=U(x _(*))−U(x _(#))   (83)which is the same form as Equation (82). U(x)=U(x) is then defined.Consequently, the location of the displacement vector U must correspondto the 0-pixels of K^(p). The previous definitions are extended bylinearity to the 1-chains of K^(p).The Force-Stress Relation

Let us consider another 3-complex K^(s) also representing the subdividedspatial support of the body B. Let us also consider a 3-cochain F and a2-cochain S such that F is the coboundary of S: $\begin{matrix}{{{\mathcal{F}:{\mathcal{C}_{3}\left( \mathcal{K}^{p} \right)}}->{\mathbb{R}}^{3}}\quad{\left. \gamma\mapsto{< \mathcal{F}} \right.,{\gamma>= < {\delta\quad\mathcal{S}}},{\gamma>= < \mathcal{S}},{{\partial\gamma} >}}} & (84)\end{matrix}$

FIGS. 48 a and 48 b present some examples of F and S for a 3-pixel ofK^(p).

Let γ_(F) be a 3-pixel of K^(s) and γ_(S) be a 2-chain over K^(s) suchthat γ_(S)=∂γ_(F). Let us assume that the 2-faces γ_(S), of γ_(F) arerelatively positively oriented with respect to the orientation of γ_(F).The definition of the coboundary leads to: $\begin{matrix}{{< \mathcal{F}},{\gamma_{F}>={\sum\limits_{\gamma_{s_{j}}}{< \mathcal{S}}}},{\gamma_{S_{j}} >}} & (85)\end{matrix}$

Again, the computational rules associated with F and S are determined.Setting a zero velocity field in Equation (71) leads to:${\underset{V}{\int\int}\rho\quad b_{i}{\mathbb{d}V}} = {\underset{S}{\int\int} - {{\sigma_{i} \cdot n}{\mathbb{d}S}}}$where σ₁=(σ_(1i), σ_(2i), σ_(3i)). To fulfill Equation (85), thefollowing relation is defined: $\begin{matrix}{{< \mathcal{F}_{i}},{\gamma_{F}>={\underset{V}{\int{\int\int}}\rho\quad b_{i}{\mathbb{d}V}\quad{and}}}} & (86) \\{{< \mathcal{S}_{i}},{\gamma_{S_{j}}>={\underset{S}{\int\int} - {{\sigma_{i} \cdot n}{\mathbb{d}S}\quad\left( {{i = 1},2,3} \right)}}}} & (87)\end{matrix}$whereF=(F ₁ , F ₂ , F ₃)^(T)andS=(S ₁ , S ₂ , S ₃)^(T)The Stress-Strain Relation

Exact global versions of Equations (66) on K^(p) and (74) on K^(s)having been presented by the means of Equations (83), (86) and (87). Inorder to complete the scheme of FIG. 44, representation of the Hooke'slaw (Equations 69 and 70) which links the local values of ε(x) and σ(x)is needed. Since Equations (69) and (70) are constitutive equations, atopological expression thereof cannot be provided. Instead, they areexpressed as linear transformations between the cochains D and S:D

S

To find this transformation, Equation (87) is recalled: < 𝒮_(i), γ_(S_(j)) >  = ∫∫_(S)−σ_(i) ⋅ n𝕕S  (i = 1, 2, 3)which links the cochain S with the strains ε using the generalizedHooke's law. Unfortunately, the strains are only known at a finitenumber of points and are then approximated over the whole domain S withan approximation function {tilde over (ε)}(x). {tilde over (ε)}(x) ischosen such that for each 1-face γ_(D) of a 2-pixel γ of K^(p), thefollowing relation is obtained: $\begin{matrix}{{{\int_{\gamma_{D}}{{\overset{\sim}{ɛ}(x)} \cdot {\mathbb{d}R}}} = {< \mathcal{D}}},{\gamma_{D} >}} & (88)\end{matrix}$where dR is an infinitesimal part of the domain represented by γ_(D). Itshould be pointed out that that only approximation of the normalcomponents of ε is needed. In fact, given:${\nabla{\overset{\sim}{U}(x)}} = {\left( {{\frac{\partial{\overset{\sim}{U}}_{1}}{\partial x_{1}}(x)},{\frac{\partial{\overset{\sim}{U}}_{2}}{\partial x_{2}}(x)},{\frac{\partial{\overset{\sim}{U}}_{3}}{\partial x_{3}}(x)}} \right)^{T} = {\left( {{{\overset{\sim}{ɛ}}_{11}(x)},{{\overset{\sim}{ɛ}}_{22}(x)},{{\overset{\sim}{ɛ}}_{33}(x)}} \right)^{T} = {\overset{\sim}{ɛ}(x)}}}$where ũ(x) is the approximated displacement vector over γ and if {tildeover (e)}(x) is chosen to satisfy Equation (88), then the vector Ũ(x) isfully determined. Then the shear components of {tilde over (ε)}(x) canbe computed by appropriately differentiating the components of Ũ(x).Using this remark and applying the generalized Hooke's law to {tildeover (ε)}(x) satisfying Equation (88), the following relation isobtained:${{\overset{\sim}{\sigma}}_{ii}(x)} = {\frac{E}{\left( {1 + v} \right)\left( {1 - {2v}} \right)}\left\lbrack {{\left( {1 - v} \right){{\overset{\sim}{ɛ}}_{ii}(x)}} + {v\left( {{{\overset{\sim}{ɛ}}_{jj}(x)} + {{\overset{\sim}{ɛ}}_{kk}(x)}} \right)}} \right\rbrack}$${{\overset{\sim}{\sigma}}_{ij}(x)} = {\frac{E}{\left( {1 + v} \right)}{{\overset{\sim}{ɛ}}_{ij}(x)}\quad\left( {i,j,{k = 1},2,3} \right)}$at all point of γ. Equation (87) is then replaced by: $\begin{matrix}{{< \mathcal{S}_{i}},{{\gamma_{S_{j}}>={\int{\int_{S}{{{- {{\overset{\sim}{\sigma}}_{i}(x)}} \cdot n}{\mathbb{d}S}}}}}\quad = {{\Lambda_{i}(ɛ)}\quad\left( {{i = 1},2,3} \right)}}} & (89)\end{matrix}$which depends on the choice of the approximation function {tilde over(ε)}(x) and of the relative position of K^(s) with respect to K^(p).Summary of the Algorithm

The algorithm used to find an expression of the internal forcesaccording to the displacements of a body is now summarized. The inputdata for this algorithm are the cochain U and the material properties ofthe body (the values of E and v).

-   -   1. Choice of the location of K_(p) with respect to K^(s).    -   2. Computation of the cochain D.    -   3. Computation of the cochain S        -   (a) Choice of the approximation function {tilde over            (ε)}(x).        -   (b) Application of Equation (89) to express S as a function            of the displacement components.    -   4. Computation of the force by applying Equation (85).        Applications        Active Contours

The above described approach is applied to a 2D active contour modelbased upon a Lagrangian evolution of the curve S: $\begin{matrix}{{\frac{\partial S}{\partial t} + {KS}} = {F_{ext}(S)}} & (90)\end{matrix}$where K is the matrix which contains the regularization forces of thecurve.

This dynamic system is discretized in time using a finite differencescheme. For a given time step Δt, the time derivative can beapproximated by:$\frac{\partial S}{\partial t} \cong \frac{S_{t + {\Delta\quad t}} - S_{t}}{\Delta\quad t}$

The curve deformation is governed by each vertex displacements comparedto their neighbors until the equilibrium between the inertia forces, theimage forces and the internal forces. Equation 90 is solved using anexplicit scheme:S _(t+Δt) =S _(t) +Δt(F _(ext) −KS _(t))   (91)

Assuming that the initial curve S₀ was in an equilibrium state and thatthe initial body forces F₀=KS₀ are constant during the deformationprocess, these forces can be added to the external forces F_(ext)leading to a modified version of Equation (91): $\begin{matrix}\begin{matrix}{S_{t + {\Delta\quad t}} = {S_{t} + {\Delta\quad{t\left( {F_{ext} + F_{0} - {KS}_{t}} \right)}}}} \\{= {S_{t} + {\Delta\quad{t\left( {F_{ext} - \left( {{KS}_{t} - {KS}_{0}} \right)} \right)}}}} \\{= {S_{t} + {\Delta\quad{t\left( {F_{ext} - {KU}} \right)}}}}\end{matrix} & (92)\end{matrix}$where U is the displacement vector of the curve S.

The image subdivision process is similar to the one presented in [57].Here, it is desired to solve Equation (92) for local U(x) located at thecenter of each pixel and known initial curve S₀ closed to the solution.Following the steps presented hereinabove, the two dimensional cubicalcomplexes K^(p) and K^(s) are first positioned. As mentioned, K^(p) isplaced in order to have its 0-pixels corresponding to the center of theimage pixels. K^(s) is positioned in such a way that its 2-pixelscoincide with the image pixels. Thus, the 2-pixels of K^(s) arerectangular and symmetrically staggered with the 1-pixels of K^(p) andeach 1-pixel of intersects orthogonally in the middle of a 1-pixel ofK^(p). Mattiussi [61] showed that this way of positioning K^(s) allowsthe use of lower order approximation polynomials without losingaccuracy. FIG. 49 shows the two complexes positions for a 5×5 image.

In order to solve Equation (92), global values F are needed over eachpixel of K^(s). Since these values are generally known, there is no needto try to express them in a topological way. In the examples, thegradient field of the bright line plausibility image obtained using aline detector proposed in [54] has been used. It was assumed that thegradient provides global values valid over the whole pixel. Thus F=ΔL isset where L is the line plausibility image, ΔL=g′_(σ) *L and g′_(σ) isthe Gaussian derivative at scale σ. An approximation function {tildeover (ε)}(x)=({tilde over (ε)}₁₁ (x), {tilde over (ε)}₂₂ (x))^(T) isalso chosen. For simplicity, it is assumed that {tilde over(ε)}(x)=∇Ũ(x) arises from a bilinear approximation:Ũ(x)=(Ũ ₁(x), Ũ ₂(x))^(T) =a+bx ₁ +cx ₂ +dx ₁ x ₂Thus:{tilde over (ε)}₁₁(x)=b+dx ₂{tilde over (ε)}₂₂(x)=c+dx ₂

Since {tilde over (ε)}₁₁ (x) and {tilde over (ε)}₂₂ (x) has to satisfyEquation (88) for all 1-faces γ_(D) of any 2-pixel γ of K^(p) as in FIG.50, the following relations hold:$\mathcal{D}_{1} = {\int_{0}^{\Delta}{{{\overset{\sim}{ɛ}\left( {x_{1},0} \right)} \cdot \overset{->}{i}}{\mathbb{d}x_{1}}}}$$\mathcal{D}_{2} = {\int_{0}^{\Delta}{{{\overset{\sim}{ɛ}\left( {\Delta,x_{2}} \right)} \cdot \overset{->}{j}}{\mathbb{d}x_{2}}}}$$\mathcal{D}_{3} = {\int_{0}^{\Delta}{{{\overset{\sim}{ɛ}\left( {x_{1},\Delta} \right)} \cdot \overset{->}{i}}{\mathbb{d}x_{1}}}}$$\mathcal{D}_{4} = {\int_{0}^{\Delta}{{{\overset{\sim}{ɛ}\left( {0,\Delta} \right)} \cdot \overset{->}{j}}{\mathbb{d}x_{2}}}}$from which the following relation is obtained: $\begin{matrix}{{\overset{\sim}{ɛ}(x)} = {{\frac{1}{\Delta}\left\lbrack {\left( {\mathcal{D}_{1} + {\frac{\mathcal{D}_{3} - \mathcal{D}_{1}}{\Delta}x_{2}}} \right),\left( {\mathcal{D}_{4} + {\frac{\mathcal{D}_{2} - \mathcal{D}_{4}}{\Delta}x_{1}}} \right)} \right\rbrack} = {\nabla{\overset{\sim}{U}(x)}}}} & (93)\end{matrix}$

From Equation 93 and the definition of the normal strains, it isstraightforward that: $\begin{matrix}{{\overset{\sim}{U}(x)} = {k + {\frac{1}{\Delta}\left( {{\mathcal{D}_{1}x_{1}} + {\mathcal{D}_{4}x_{2}}} \right)} + {\frac{1}{\Delta^{2}}\left( {\mathcal{D}_{3} - \mathcal{D}_{1} + \mathcal{D}_{2} - \mathcal{D}_{4}} \right)x_{1}x_{2}}}} & (94)\end{matrix}$where k=Ũ(0). Equations (83) and (94) lead to: $\begin{matrix}\begin{matrix}{{\overset{\sim}{U}(x)} = {\overset{\sim}{U}\left( {x_{1},x_{2}} \right)}} \\{= {{U\left( {0,0} \right)} + {\frac{1}{\Delta}\left( {{U\left( {\Delta,0} \right)} - {U\left( {0,0} \right)}} \right)x_{1}} + {\frac{1}{\Delta}\left( {{U\left( {0,\Delta} \right)} -} \right.}}} \\{{\left. {U\left( {0,0} \right)} \right)x_{2}} + {\frac{1}{\Delta^{2}}\left( {{U\left( {0,0} \right)} + {U\left( {\Delta,\Delta} \right)} - {U\left( {0,\Delta} \right)} -} \right.}} \\{\left. {U\left( {\Delta,0} \right)} \right)x_{1}x_{2}}\end{matrix} & (95)\end{matrix}$from which the values of {tilde over (σ)}_(i) (x)=({tilde over(σ)}_(1i)′(x), {tilde over (σ)}_(2i) (x))^(T) can be deduced.

The last step is the computation of the internal forces F for each2-pixel of K^(s). With K^(p) and K^(s) positioned as mentioned before,each 2-pixel γ_(F) of K^(s) intersects four 2-pixels γ_(A), γ_(B), γ_(C)and γ_(D) of K^(p). That is, four approximation functions {tilde over(σ)}_(i) ^(A), {tilde over (σ)}_(i) ^(B), {tilde over (σ)}_(i) ^(C) and{tilde over (σ)}_(i) ^(D) corresponding to the four intersecting2-pixels of K^(p) (see FIG. 51) have to be considered.

The value of the cochain S are found over the four 1-face of γ by theappropriate integration:$\mathcal{S}_{i}^{1} = {{\int_{0}^{\frac{\Delta}{2}}{{{{\overset{\sim}{\sigma}}_{i}^{A}\left( {\frac{\Delta}{2},x_{2}} \right)} \cdot \overset{->}{i}}{\mathbb{d}x_{2}}}} + {\int_{- \frac{\Delta}{2}}^{0}{{{{\overset{\sim}{\sigma}}_{i}^{B}\left( {\frac{\Delta}{2},x_{2}} \right)} \cdot \overset{->}{i}}{\mathbb{d}x_{2}}}}}$$\mathcal{S}_{i}^{2} = {{\int_{0}^{\frac{\Delta}{2}}{{{{\overset{\sim}{\sigma}}_{i}^{B}\left( {x_{1},{- \frac{\Delta}{2}}} \right)} \cdot {- \overset{\rightarrow}{j}}}{\mathbb{d}x_{1}}}} + {\int_{- \frac{\Delta}{2}}^{0}{{{{\overset{\sim}{\sigma}}_{i}^{C}\left( {x_{1},{- \frac{\Delta}{2}}} \right)} \cdot {- \overset{\rightarrow}{j}}}{\mathbb{d}x_{1}}}}}$$\mathcal{S}_{i}^{3} = {{\int_{0}^{\frac{\Delta}{2}}{{{{\overset{\sim}{\sigma}}_{i}^{A}\left( {x_{1},\frac{\Delta}{2}} \right)} \cdot \overset{->}{j}}{\mathbb{d}x_{1}}}} + {\int_{- \frac{\Delta}{2}}^{0}{{{{\overset{\sim}{\sigma}}_{i}^{D}\left( {x_{1},\frac{\Delta}{2}} \right)} \cdot \overset{->}{j}}{\mathbb{d}x_{1}}}}}$$\mathcal{S}_{i}^{4} = {{\int_{0}^{\frac{\Delta}{2}}{{{{\overset{\sim}{\sigma}}_{i}^{D}\left( {{- \frac{\Delta}{2}},x_{2}} \right)} \cdot {- \overset{->}{i}}}{\mathbb{d}x_{2}}}} + {\int_{- \frac{\Delta}{2}}^{0}{{{{\overset{\sim}{\sigma}}_{i}^{C}\left( {{- \frac{\Delta}{2}},x_{2}} \right)} \cdot {- \overset{->}{i}}}{\mathbb{d}x_{2}}}}}$Equation (85) leads to:<F _(i), γ_(F) >=S _(i) ¹ +S _(i) ² +S _(i) ³ +S _(i) ⁴   (96)

By substituting Equations (96) and (95) in Equation (96), the internalforces F can be expressed as a function of the displacement U. As anexample, the values of F are represented for the 2-pixel γ_(F) of FIG.50 with Δ=1:F₁ = C[(3 − 4v)u_(−1, 1) + (2 − 8v)u_(0, 1) + (3 − 4v)u_(1, 1) + (10 + 8v)u_(−1, 0) + (−36 + 48v)u_(0, 0) + (10 − 8v)u_(1, 0) + (3 − 4v)u_(−1, −1) + (2 − 8v)u  0, −1 + (3 − 4v)u_(1, −1) − 2υ_(−1, 1) + 2υ_(1, 1) + 2υ_(−1, −1) − 2υ_(1, −1)]F₂ = C[(3 − 4v)u_(−1, 1) + (10 − 8v)u_(0, 1) + (3 − 4v)u_(1, 1) + (2 − 8v)u_(−1, 0) + (−36 + 48v)u_(0, 0) + (2 − 8v)u_(1, 0) + (3 − 4v)u_(−1, −1) + (10 − 8v)u_(0, −1) + (3 − 4v)u_(1, −1) − 2υ_(−1, 1) + 2υ_(1, 1) + 2υ_(−1, −1) − 2υ_(1, −1)]${{{where}\text{:}\quad C} = {{\frac{E}{16\left( {1 + v} \right)\left( {1 - {2v}} \right)}\quad{and}\quad\mathcal{U}} = \begin{bmatrix}u \\v\end{bmatrix}}}\quad$

Equation (97) induces a linear relationship between a pixel and itsneighbors. This relation is used to build the stiffness matrix ofEquation (91). If U_(x) ₁ (i=1, 2) is considered as the displacementvector for the component x_(i) then:F _(i)(x)=U _(x) ₁ *N _(x) ₁ ^(i) +U _(x) ₂ *N _(x) ₂ ^(i)where$N_{x_{1}}^{i} = {\frac{E}{16\left( {1 + v} \right)\left( {1 - {2v}} \right)}\begin{bmatrix}{3 - {4v}} & {2 - {8v}} & {3 - {4v}} \\{10 - {8v}} & {{- 36} + {48v}} & {10 - {8v}} \\{3 - {4v}} & {2 - {8v}} & {3 - {4v}}\end{bmatrix}}$$N_{x_{2}}^{i} = {\frac{E}{16\left( {1 + v} \right)\left( {1 - {2v}} \right)}\begin{bmatrix}{- 2} & 0 & {2} \\{0} & 0 & {0} \\2 & 0 & {- 2}\end{bmatrix}}$$N_{x_{3}}^{i} = {\frac{E}{16\left( {1 + v} \right)\left( {1 - {2v}} \right)}\begin{bmatrix}{- 2} & 0 & {2} \\{0} & 0 & {0} \\2 & 0 & {- 2}\end{bmatrix}}$$N_{x_{4}}^{i} = {\frac{E}{16\left( {1 + v} \right)\left( {1 - {2v}} \right)}\begin{bmatrix}{3 - {4v}} & {10 - {8v}} & {3 - {4v}} \\{2 - {8v}} & {{- 36} + {48v}} & {2 - {8v}} \\{3 - {4v}} & {10 - {8v}} & {3 - {4v}}\end{bmatrix}}$

The pairs (N_(x) _(i) ^(l), N_(x) ₂ ^(l)), (i=1, 2) will be referred toas the stiffness kernels.

Computation of the Displacement Vector

The assumption made when calculating the displacement vector will now beexplained. Let v be a vertex of subdivided curve S and v′ be itscorresponding vertex in the deformed curve S′. Let us denote by U[v] theentry in the displacement vector U corresponding to v. Let us supposethat the displacement is constant in each direction everywhere that isU[v]=(k₁, k₂)^(T) with k₁, k₂ ε R for all v in S. Since the sum of allentries of either N_(x) ₁ ¹, N_(x) ₂ ¹, N_(x) ₁ ² or N_(x) ₂ ² is zero,it follows that F₁[v]=F₂[v]=0 for all v in S which means that there isno internal force induced. The computation of the internal forces withthe stiffness kernels is then invariant with respect to translation.

Let v₁, v₂, v₃, v₄ and v₅ be five adjacent vertices of S and v′₁, v′₂,v′₃, v′₄, v′₅ be their corresponding vertices in S′ (see FIG. 52).

Let v_(i,x) _(j) be the x_(j) coordinate of the spatial representationof the vertex v_(i). Then, the following relation is obtained:U _(x) _(i) =( . . . , v′ _(1,x) _(i) −v _(1,x) _(i) , v′ _(2,x) _(i) −v_(2,x) _(i) , v′ _(3,x) _(i) −v _(3,x) _(i) , v′ _(4,x) _(i) −v _(4,x)_(i) , v′ _(5,x) _(i) −v _(5,x) _(i) , . . . )^(T) (i=1

The translation invariance property leads to: $\begin{matrix}{{F_{i}(x)} = {{U_{x_{1}}*N_{x_{1}}^{i}} + {U_{x_{2}}*N_{x_{2}}^{i}}}} \\{= {{\left\lbrack {U_{x_{1}} - \left\lbrack {\upsilon_{2,x_{1}}^{\prime} - \upsilon_{2,x_{1}}} \right\rbrack} \right\rbrack*N_{x_{1}}^{i}} + {\left\lbrack {U_{x_{2}} - \left\lbrack {\upsilon_{2,x_{2}}^{\prime} - \upsilon_{2,x_{2}}} \right\rbrack} \right\rbrack*N_{x_{2}}^{i}}}}\end{matrix}$where [v′_(2,x) ₁ −v_(2,x) ₁ ] stands for a matrix whose all entriesequal v′_(2,x) ₁ −v_(2,x) ₁ . The displacement component used to computethe internal force F_(i) at vertex v₃ is then:U_(x_(k))[v₃] = (v_(3, x_(k))^(′) − v_(3, x_(k))) − (v_(2, x_(k))^(′) − v_(2, x_(k)))   = (v_(3, x_(k))^(′) − v_(2, x_(k))^(′)) − (v_(3, x_(k)) − v_(2, x_(k)))which is the relative displacement of the vertex v₃ with respect to v₂.However, nothing prevents computing the relative displacement of v₃ withrespect to v₃. In this case, the x_(k) displacement component used tocompute the internal force F_(i) at vertex v₃ would be: $\begin{matrix}{{U_{x_{k}}\left\lbrack v_{3} \right\rbrack} = {\left( {v_{3,x_{k}}^{\prime} - v_{3,x_{k}}} \right) - \left( {v_{4,x_{k}}^{\prime} - v_{4,x_{k}}} \right)}} \\{= {\left( {v_{3,x_{k}}^{\prime} - v_{4,x_{k}}^{\prime}} \right) - \left( {v_{3,x_{k}} - v_{4,x_{k}}} \right)}}\end{matrix}$

In order to take into account these facts, the average value of therelative displacements for all adjacent vertices to v₂ is used. Thus:$\begin{matrix}{{U_{x_{k}}\left\lbrack v_{2} \right\rbrack} = {\frac{1}{2}\begin{bmatrix}{\left( {v_{3,x_{k}}^{\prime} - v_{2,x_{k}}^{\prime}} \right) - \left( {v_{3,x_{k}} - v_{2,x_{k}}} \right) +} \\{\left( {v_{3,x_{k}}^{\prime} - v_{4,x_{k}}^{\prime}} \right) - \left( {v_{3,x_{k}} - v_{4,x_{k}}} \right)}\end{bmatrix}}} & (97)\end{matrix}$

Reorganizing the terms of Equation (97), for an arbitrary vertex v_(i)of S, the following relation is obtained: $\begin{matrix}\begin{matrix}{{U_{x_{k}}\left\lbrack v_{i} \right\rbrack} = {\left\lbrack \frac{v_{{i + 1},x_{k}} - {2v_{i,x_{k}}} + v_{{i - 1},x_{k}}}{2} \right\rbrack -}} \\{\left\lbrack \frac{v_{{i + 1},x_{k}}^{\prime} - {2v_{i,x_{k}}^{\prime}} + v_{{i - 1},x_{k}}^{\prime}}{2} \right\rbrack\quad{or}} \\{{U\left\lbrack v_{i} \right\rbrack} = {\left\lbrack \frac{v_{i + 1} - {2v_{i}} + v_{i - 1}}{2} \right\rbrack - \left\lbrack \frac{v_{i + 1}^{\prime} - {2v_{i}^{\prime}} + v_{i - 1}^{\prime}}{2} \right\rbrack}} \\{= {\frac{1}{2}\left\lbrack {\frac{\partial^{2}S}{\partial v^{2}} - \frac{\partial^{2}S^{\prime}}{\partial v^{2}}} \right\rbrack}}\end{matrix} & (98)\end{matrix}$if we assume a finite difference approximation of the second derivativeof S and S′ with respect to their vertices.

Let us finally notice that the second derivative of the curve S inEquation (98) can be computed using Gaussian derivatives:$\frac{\partial^{2}S}{\partial v^{2}} = {g_{\sigma}^{''}*S}$where σ controls the degree of smoothing. Such a computation of thesecond derivative of S allows to obtain smooth results by simulating asmoother target curve.Experimental ResultsActive Contours

The approach proposed herein has been experimented on real and syntheticimages in the context of high-resolution images of road databases. Foreach image, the results have been compared with another method.

FIG. 55 a presents the results for the physics-based method (PB)according to the present invention for an aerial image while FIG. 55 bshows the results for the finite element (FEM) method (α and β unknown).A material similar to rubber (see the Table in FIG. 53) with E=150 andv=0.45 has been simulated. In both images, the image force is thegradient (σ=1.5) of the bright line plausibility image obtained using aline detector proposed in [54] with the line detection scale set to 0.8.FIGS. 54 a and 54 b show respectively the initial curve S₀ and thebright line plausibility image.

FIG. 56 shows a SAR image in which the initial curves are drawn inwhite. The line plausibility image (σ=1.5 and line detection scale=0.8)of both bright and dark lines is shown in FIG. 57. FIGS. 58 and 59presents respectively the results obtained with PB (E=150 and v=0.45)and FEM (α and β unknown) methods. One can notice that the PB curves arecloser to the shore than the FEM especially in region of high curvature.

FIG. 60 shows some initialization for the first band of a Landsat 7image. The bright line plausibility image (σ=1.5 and line detectionscale=0.8) is shown in FIG. 61. FIGS. 62 and 63 present the resultsobtained with the PB and FEM methods respectively.

The fact that the deformations obtained using the model according tothis illustrative embodiment of the present invention have a physicalinterpretation has been discussed. The fact that the objects modeledusing the PB method have their own physical properties and the abilityto recover their original shape when the external forces applied on themare removed has also been discussed. To illustrate this fact, FIGS. 64 aand 64 b show some initialization and the corrected curve for asynthetic image. FIGS. 65 a through FIGS. 65 d show the evolution ofthis corrected curve when the external forces are removed. FIG. 65 dpresents both the final curve (in black) and the initial curve (inwhite). One can clearly notice that the curve has recovered its originalshape but has also experienced a spatial shift.

Conclusion

A new model for active contours was presented. The proposed approachdecompose the image using an image model based on algebraic topology.This model uses generic mathematical tools which can be applied to solveother problems such as linear and non-linear diffusion and optical flow[57]. Moreover, the model works with either 2D or 3D images and caneasily be extended to active surfaces and active volumes.

The approach presents the following differences with the othermethods: 1) Both global and local quantities are used; 2) The model isbased upon basic laws of physics. This allows us to give a physicalexplanation to the deformation steps; 3) The curves and surfaces havephysical behaviors such as the ability to recover their original shapeonce the applied forces are removed; 4) Approximation are made only whenthe constitutive equation is involved.

Although the present invention has been described hereinabove by way ofnon-restrictive illustrative embodiments thereof, it can be modified atwill, within the scope of the appended claims, without departing fromthe spirit and nature of the subject invention.

REFERENCES

-   [1] M. Allili and D. Ziou, Extraction of topological properties of    images via cubical homology, Technical Report, 2000.-   [2] M. Allili and D. Ziou, Topological Feature Extraction in Binary    Images, IEEE ISSPA, Malysia, August 2001.-   [3] M. Allili, K. Mischaikow, and A. Tannenbaum, Cubical Homology    and the Topological Classification of 2D and 3D Imagery, Int. Conf.    Image Processing, Greece, 2001.-   [4] M. F. Auclair-Fortier, P. Poulin, D. Ziou, and M. Allili,    Computational Algebraic Topology for Resolution of the Poisson    Equation: Application to Computer Vision, Technical Report, 2001.-   [5] R. Egli and N. F. Stewart, A framework for system specification    using chains on cell complexes, Computer Aided Design 32, 447-459,    2000.-   [6] P. J. Giblin, Graphs, Surfaces and Homology, Chapman and Hall,    London, 1977.-   [7] T. Y. Kong and A. Rosenfeld, Digital Topology: Introduction and    Survey, CVGIP 48, 357-393, 1989.-   [8] V. A. Kovalvesky, Finite Topology to Image Analysis, Comp.    Vision, and Image Processing 46, 141-161, 1989.-   [9] W. S. Massey, A Basic Course in Algebraic Topology,    Springer-Verlag, New York, 1991.-   [10] J. R. Munkres, Elements of Algebraic Topology, Addison-Wesley,    1984.-   [11] R. S. Palmer and V. Shapiro, Chain Models of Physical Behavior    for Engineering Analysis and Design, Research in Engineering Design    5, 161-184, 1993.-   [12] P. Poulin, D. Ziou, and M. F. Auclair-Fortier, Computational    Algebraic Topology for Deformable Objects, Technical Report, 2001.-   [13] A. S. Schwarz, Topology for Physicists, Springer-Verlag, 1994.-   [14] E. Tonti, On the Formal Structure of Physical Theories,    Technical Report Istituto Di Mathematica Del Politechnico Di Milano,    Milan, 1975.-   [15] D. Ziou and M. Allili, Computation of the Euler Number in    Binary Images without Skeletonization via Cubical Complex, Technical    Report, 2001.-   [16] M. Allili and D. Ziou. Extraction of Topological Properties of    Images via Cubical Homology. Technical Report CDSNS 2000-365,    Georgia Institute of Technology, June 2000.-   [17] L. Alvarez, R. Deriche, and F. Santana. Recursivity and PDE's    in Image Processing. In Proceedings 15th International Conference on    Pattern Recognition, volume I, pages 242-248, September 2000.-   [18] J. L. Barron, D. J. Fleet, and S. S. Beauchemin. Systems and    Experiment Performance of Optical Flow Techniques. International    Journal of Computer Vision, 12(1):43-77, 1994.-   [19] S. S. Beauchemin and J. L. Barron. The Computation of Optical    Flow. ACM Computing Survey, 27(3), 1995.-   [20] A. J. Chapman. Fundamentals of Heat Transfer. Macmillan    Publishing Company, 1987.-   [21] A. K. Chhabra and T. A. Grogan. On Poisson Solvers and    Semi-Direct Methods for Computing Area Based Optical Flow. IEEE    Transactions on Pattern Analysis and Machine Intelligence,    16:1133-1138, 1994.-   [22] H. Chidiac and D. Ziou. Formation d'images optiques. Technical    Report 226, Département de mathématiques et d'informatique,    Universitë de Sherbrooke, November 1998.-   [23] L. D. Cohen and I. Cohen. Finite Element Methods for Active    Contour Models and Balloons for 2D and 3D Images. IEEE Transactions    on Pattern Analysis and Machine Intelligence, 15, Nov. 1993.-   [24] L. D. Cohen. On active contour models and balloons. Computer    Vision, Graphics and Image Processing, 53(2):211-218, March 1991.-   [25] R. Deriche and O. Faugeras. Les EDP en traitement des images et    vision par ordinateur. Technical Report 2697, INRIA, November 1995.-   [26] C. H. Edwards and D. E. Penney. Calculus with Analytic    Geometry. Prentice Hall, 1998.-   [27] H. U. Fuchs. The Dynamics of Heat. Springer-Verlag, 1996.-   [28] D. Halliday, R. Resnick, and J. Walker. Fundamentals of    Physics. John Willey and Sons, 1997.-   [29] B. K. P. Horn and B. Schunck. Determining Optical Flow.    Artificial Intelligence, 17:185-203, 1981.-   [30] J. Kervorkian. Partial Differential Equations: Analytical    Solution Techniques, chapter 2, pages 48-116. Mathematics Series.    Chapman and Hall, 1990.-   [31] S. H. Lai and B. C. Vermuri. An O(n) Iterative Solution to the    Poisson Equation in Low-level Vision Problems. Technical Report    TR93-035, University of Florida, Computer and Information Sciences    Department, 1993.-   [32] J. Li and A. O. Hero. A Spectral Method for Solving Elliptic    Equations for Surface Reconstruction and 3D Active Contours.    Proceedings of IEEE International Conference on Image Processing,    Thessaloniki, Greece, October 2001.-   [33] G. T. Mase and G. E. Mase. Continuum Mechanics for Engineers.    CRC Press, 1999.-   [34] C. Mattiussi. An Analysis of Finite Volume, Finite Element, and    Finite Difference Methods Using Some Concepts from Algebraic    Topology. Journal of Computational Physics, 133:289-309, 1997.-   [35] J. Monteil and A. Beghdadi. A New Interpretation and    Improvement of the Nonlinear Anisotropic Diffusion for Image    Enhancement. IEEE Transactions on Pattern Analysis and Machine    Intelligence, 21(9):940-946, September 1999.-   [37] S. V. Patankar. Numerical Heat Transfer and Fluid Flow.    Computational Methods in Mechanics and Thermal Sciences. McGraw-Hill    Book Company, 1980.-   [38] P. Perona and J. Malik. Scale-Space and Edge Detection Using    Anisotropic Diffusion. IEEE Transactions on Pattern Analysis and    Machine Intelligence, 12(7):629-639, July 1990.-   [39] T. Symchony, R. Chellappa, and M. Shao. Direct Analytical    Methods for Solving Poisson Equations in Computer Vision Problems.    IEEE Transactions on Pattern Analysis and Machine Intelligence,    12(5):435-446, May 1990.-   [40] D. Terzopoulos. Image Analysis Using Multigrid Relaxation    Methods. IEEE Transactions on Pattern Analysis and Machine    Intelligence, 8:129-129, 1986.-   [41] G. B. Thomas and R. L. Finney. Calculus and Analytic Geometry.    Addison-Wesley Publishing Company, 1988.-   [42] E. Tonti. On the Formal Structure of Physical Theories.    Technical report, Istituto Di Mathematica Del Polithechnico Di    Milano, 1975.-   [43] J. Weickert. A Review of Nonlinear Diffusion Filtering. Lecture    Notes in Computer Science, 1252:3-28, 1997.-   [44] E. Zauderer. Partial Differential Equations of Applied    Mathematics. Pure and Applied Mathematics. John Willey and Sons, New    York, US, 1983.-   [45] M. Allili and D. Ziou. Extraction of Topological Properties of    Images via Cubical Homology. Technical Report CDSNS 2000-365,    Georgia Institute of Technology, June 2000.-   [46] M.-F. Auclair-Fortier, D. Ziou, C. Armenakis, and S. Wang.    Automated Correction and Updating of Road Databases from    High-Resolution Imagery. Canadian Journal of Remote Sensing,    27:76-89, 2001.-   [48] K. J. Bathes. Finite element procedures. Prentice Hall, 1996.-   [49] A. P. Boresi. Elasticity in Engineering Mechanics. Prentice    Hall, 1965.-   [50] A. P. Boresi, R. J. Schmidt, and O. M. Sidebottom. Advanced    Mechanics of Materials, Fifth edition. John Willey and Sons, 1993.-   [51] M. Bro-Nielsen. Surgery simulation using fast finite elements.    In Proc. Visualization in Biomedical Computing (VBC'96), volume    1131, pages 529-534, Hamburg, Germany, September 1996. Springer    Lecture Notes in Computer Science.-   [52] E. F. Byars and R. D. Snyder. Mechanics of Deformable Bodies.    Intext Educational Publishers, 1975.-   [53] S. Cotin and N. Ayache. A Hybrid Elastic Model Allowing    Real-Time Cutting, Deformations and Force-Feedback for Surgery    Training and Simulation. In Proc. of CAS99, pages 70-81, May 1999.-   [54] F. Deschénes, D. Ziou, and M.-F. Auclair-Fortier. Detection of    Lines, Line Junctions and Line Terminations. Technical Report 259,    Département de mathématiques et d'informatique, Université de    Sherbrooke, 2000. Submitted in ISPRS Journal of Photogrammetry and    Remote Sensing, 2000.-   [55] C. H. Edwards and D. E. Penney. Calculus with Analytic    Geometry. Prentice Hall, 1998.-   [56] K. M. Entwistle. Basic Principles of the Finite Element Method.    IOM Communications Ltd, 1999.-   [57] S. F. F. Gibson and B. Mirtich. A Survey of Deformable Modeling    in Computer Graphics. Technical report, Mitsubishi Electric Research    Laboratory, 1997.-   [58] R. C. Juvinall. Engineering considerations of Stress, Strain    and Strength. McGraw-Hill, 1967.-   [59] M. Kass, A. Witkin, and D. Terzopoulos. Snakes: Active Contour    Models. The International Journal of Computer Vision, 1(4):321-331,    1988.-   [60] G. T. Mase and G. E. Mase. Continuum Mechanics for Engineers.    CRC Press, 1999.-   [61] C. Mattiussi. An Analysis of Finite Volume, Finite Element, and    Finite Difference Methods Using Some Concepts from Algebraic    Topology. Journal of Computational Physics, 133:289-309, 1997.-   [62] J. Montagnat and H. Delingette. Volumetric Medical Images    Segmentation using Shape Constrained Deformable Models. Lecture    Notes in Computer Science, 1205:13-22, 1997.-   [63] J. Montagnat, H. Delingette, N. Scapel, and N. Ayache.    Representation, shape, topology and evolution of deformable    surfaces. Application to 3D medical imaging segmentation. Technical    Report 3954, INRIA, 2000.-   [64] K. W. Morton and D. F. Mayers. Numerical Solution of Partial    Differential Equations. Cambridge University Press, 1994.-   [65] B. V. Muvdi and J. W. McNabb. Engineering Mechanics of    Materials. Macmillan Publishing Co, 1980.-   [66] N. O. Myklestad. Statics of deformable bodies. MacMillan    Company, 1966.-   [67] R. S. Palmer and V. Shapiro. Chain Models of Physical Behavior    for Engineering Analysis and Design. Technical Report TR93-1375,    Cornell University, Computer Science Department, August 1993.-   [68] N. Peterfreund. Robust Tracking of Position and Velocity With    Kalman Snakes. IEEE Transactions on Pattern Analysis and Machine    Intelligence, 21(6):564-569, 1999.-   [69] W. D. Pilkey and W. Wunderlich. Mechanics of structures:    variational and computational methods. CRC Press, 1994.-   [70] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery.    Numerical Recipes in C. Cambridge University Press, 1992.-   [71] G. B. Thomas and R. L. Finney. Calculus and Analytic Geometry.    Addison-Wesley Publishing Company, 1988.-   [72] E. Tonti. On the Formal Structure of Physical Theories.    Technical report, Istituto Di Mathematica Del Polithechnico Di    Milano, 1975.

1. A computational image model, comprising: an image support including astructure of n-pixels comprising pixel faces; quantities related toimage features; and an algebraic structure relating the quantities tothe n-pixels and/or pixel faces, the algebraic structure comprisingalgebraic operations defining a relation between the quantities.
 2. Acomputational image model as defined in claim 1, wherein each n-pixel isdefined as a geometrical structure comprising vertices, edges, faces anda volume, and wherein each n-pixel comprises: a first pixel dimensionn=0 including the vertices of the n-pixel; a second pixel dimension n=1including the edges of the n-pixel; a third pixel dimension n=2including the faces of the n-pixel; a fourth pixel dimension n=3including the volume of the n-pixel; and a n^(th) pixel dimension nincluding the hypervolume of the n-pixel.
 3. A computational image modelas defined in claim 1, wherein the geometrical structure is selectedfrom the group consisting of: a cube, a triangle, a hexagone and apentagons.
 4. A computational image model as defined in claim 1, whereinthe quantities related to image features are selected from the groupconsisting of: scalar quantities, vectors, tensors and matrices.
 5. Acomputational image model as defined in claim 1, wherein the algebraicoperations comprise problem-independent operations.
 6. A computationalimage model as defined in claim 1, wherein the algebraic operationscomprise problem-dependent operations.
 7. A computational image model asdefined in claim 1, wherein the structure of n-pixels comprises pairs ofdisjoint n-pixels.
 8. A computational image model as defined in claim 1,wherein the structure of n-pixels comprises pairs of n-pixelsintersecting through a common i-pixel, where i<n.
 9. A computationalimage model as defined in claim 1, wherein each n-pixel is translatedalgebraically into a q-pixel, wherein q ε {1, 2, . . . , n}.
 10. Acomputational image model as defined in claim 9, wherein each q-pixelincludes (q−1)-faces, (q−2)-faces, . . . , (q-q)-faces.
 11. Acomputational image model as defined in claim 9, wherein the imagesupport comprises a geometrical complex, which is a collection ofq-pixels.
 12. A computational image model as defined in claim 10,wherein the image support comprises a geometrical complex, which is acollection of q-pixels, and wherein: every face of a q-pixel in thegeometrical complex is also located in the geometrical complex; and anypair of two q-pixels of the geometrical complex have an intersectionwhich is either empty or constituted by a common face of both q-pixelsof the pair.
 13. A computational image model as defined in claim 11,comprising a plurality of image supports forming the geometricalcomplex.
 14. A computational image model as defined in claim 11, whereinthe geometrical complex is expressed in algebraic form as a q-chain,which is a linear combination of all the q-pixels of the geometricalcomplex.
 15. A computational image model as defined in claim 9, whereinthe geometrical complex comprises q-cochains, which are relationsassociating quantities related to image features to the q-pixels and/orfaces of said q-pixels.
 16. A computational image model as defined inclaim 15, wherein the quantities related to image features andassociated to the q-pixels and/or faces of said q-pixels are globalquantities associated to all the q-pixels.
 17. A computational imagemodel as defined in claim 15, wherein the quantities related to imagefeatures and associated to the q-pixels and/or faces of said q-pixelsare local quantities each associated to one q-pixel and/or faces of saidone q-pixel.
 18. A computational image model as defined in claim 16,comprising (q≧1)-cochains to represent the local quantities.
 19. Acomputational image model as defined in claim 17, comprising 0-cochainto represent the global quantities.
 20. A computational image model asdefined in claim 17, wherein the algebraic operations comprise acoboundary operation giving a relationship between the q-cochains.
 21. Acomputational image model as defined in claim 9, wherein: the imagesupport comprises a plurality of geometrical complexes, each being acollection of q-pixels; and the algebraic operations comprise a codualoperation establishing a link between q-cochains that belong todifferent geometrical complexes.
 22. A method of computationallymodelling an image, comprising: producing an image support including astructure of n-pixels comprising pixel faces; defining quantitiesrelated to image features; and relating the quantities to the n-pixelsand/or pixel faces through an algebraic structure, and relating thequantities to each other through algebraic operations.
 23. A method ofcomputationally modelling an image as defined in claim 22, whereinrelating the quantities to the n-pixels and/or pixel faces through analgebraic structure comprises translating each n-pixel algebraicallyinto a q-pixel, wherein q ε {1, 2, . . . , n}, wherein each q-pixelincludes (q−1)-faces, (q−2)-faces, . . . , (q-q)-faces.
 24. A method ofcomputationally modelling an image as defined in claim 22, whereinproducing an image support comprises forming a geometrical complex,which is a collection of q-pixels, and wherein: every face of a q-pixelin the geometrical complex is also located in the geometrical complex;and any pair of two q-pixels of the geometrical complex have anintersection which is either empty or constituted by a common face ofboth q-pixels of the pair.
 25. A method of computationally modelling animage as defined in claim 24, wherein producing an image supportcomprises forming a plurality of image supports forming the geometricalcomplex.
 26. A method of computationally modelling an image as definedin claim 24, wherein relating the quantities to the n-pixels and/orpixel faces through an algebraic structure comprises expressing thegeometrical complex in algebraic form as a q-chain, which is a linearcombination of all the q-pixels of the geometrical complex.
 27. A methodof computationally modelling an image as defined in claim 24, whereinrelating the quantities to the n-pixels and/or pixel faces through analgebraic structure comprises forming, in the geometrical complex,q-cochains which are relations associating quantities related to imagefeatures to the q-pixels and/or faces of said q-pixels.
 28. A method ofcomputationally modelling an image as defined in claim 22, whereindefining quantities related to image features comprises defining globalquantities associated to all the q-pixels.
 29. A method ofcomputationally modelling an image as defined in claim 22, whereindefining quantities related to image features comprises defining localquantities associated to one q-pixel and/or faces of said one q-pixel.30. A method of computationally modelling an image as defined in claim27, wherein relating the quantities to each other through algebraicoperations comprise producing a coboundary operator giving arelationship between q-cochains.
 31. A method of computationallymodelling an image as defined in claim 27, wherein: producing an imagesupport comprises forming a plurality of geometrical complexes, eachbeing a collection of q-pixels; and relating the quantities to eachother through algebraic operations comprises producing a codualoperation establishing a link between cochains that belong to differentgeometrical complexes.
 32. An image modelling method as defined in claim27, wherein relating the quantities to the n-pixels and/or pixel facesthrough an algebraic structure comprises expressing a global quantityassociated with all q-pixels through a q-cochain such that, for twoadjacent q-pixels c_(q) ¹ and c_(q) ², the q-cochain F_(q) satisfies therelation F_(q)(λ₁c_(q) ¹+λ₂c_(q) ²)=λ₁F_(q)(c_(q) ¹)+λ₂F_(q)(c_(q) ²),where λ1 and λ2 are integers.
 33. An image modelling method as definedin claim 22, wherein: relating the quantities to the n-pixels and/orpixel faces through an algebraic structure comprises translating eachn-pixel algebraically into a q-pixel, wherein q ε {1, 2, . . . , n},wherein each q-pixel includes (q−1)-faces, (q−2)-faces, . . . ,(q-q)-faces; producing an image support comprises forming geometricalcomplexes, each being a collection of q-pixels; relating the quantitiesto the n-pixels and/or pixel faces through an algebraic structurecomprises: expressing each geometrical complex in algebraic form as aq-chain, which is a linear combination of all the q-pixels of thegeometrical complex; forming, in the geometrical complexes, q-cochainswhich are relations associating quantities related to image features tothe q-pixels and/or faces of said q-pixels; relating the quantities toeach other through algebraic operations comprises: producing acoboundary operator giving a relationship between the q-cochains; andproducing a codual operation establishing a link between q-cochains thatbelong to different geometrical complexes.
 34. A computational frameworkfor solving a problem using an image computationally modelled by meansof the method of claim 33, comprising: identifying basic laws associatedto the problem; from the identified basic laws, defining quantitiesrelated to the problem; associating the quantities to respectiveq-cochains; associating the basic laws related to the problem torespective coboundary and codual operations; and resolving the resultingalgebraic system.
 35. A computational framework as defined in claim 34,wherein forming geometrical complexes comprises forming first and secondgeometrical complexes.
 36. A computational framework as defined in claim35, wherein identifying basic laws associated to the problem comprisessupporting one basic law through the first geometrical complex.
 37. Acomputational framework as defined in claim 36, wherein the problem tobe solved is a 2D global differential equation for heat flow in ahomogeneous medium, and wherein said one basic law is a heat flow law.38. A computational framework as defined in claim 37, whereinassociating the quantities to respective q-cochains comprisesrepresenting a global quantity of temperature through a 0-cochain, andassociating the heat flow law through a 1-cochain.
 39. A computationalframework as defined in claim 35, wherein identifying basic lawsassociated to the problem comprises supporting one basic law through thesecond geometrical complex.
 40. A computational framework as defined inclaim 39, wherein the problem to be solved is a 2D global differentialequation for heat flow in a homogeneous medium, and wherein said onebasic law is a heat source law.
 41. A computational framework as definedin claim 36, wherein identifying basic laws associated to the problemcomprises supporting a second basic law through the second geometricalcomplex, and wherein associating the basic laws related to the problemto respective coboundary and codual operations comprises representing aconstitutive law linking basic laws from the first and secondgeometrical complexes by a codual operation.
 42. An image modellingmethod as defined in claim 22, wherein: relating the quantities to then-pixels and/or pixel faces through an algebraic structure comprisestranslating each n-pixel algebraically into a q-pixel, wherein q ε {1,2, . . . , n}, wherein each q-pixel includes (q−1)-faces, (q−2)-faces, .. . , (q-q)-faces; producing an image support comprises forming ageometrical complex, which is a collection of q-pixels; relating thequantities to the n-pixels and/or pixel faces through an algebraicstructure comprises: expressing the geometrical complex in algebraicform as a q-chain, which is a linear combination of all the q-pixels ofthe geometrical complex; forming, in the geometrical complex, q-cochainswhich are relations associating quantities related to image features tothe q-pixels and/or faces of said q-pixels; relating the quantities toeach other through algebraic operations comprises: producing coboundaryoperations giving a relationship between the q-cochains.
 43. Acomputational framework for solving a problem using an imagecomputationally modelled by means of the method of claim 42, comprising:identifying basic laws associated to the problem; from the identifiedbasic laws, defining quantities related to the problem; associating thequantities to respective q-cochains; associating the basic laws relatedto the problem to respective coboundary operations; and resolving theresulting algebraic system.
 44. A computational framework for solving aheat transfer problem, comprising: producing an image support includinga structure of n-pixels, the image support comprising: q-pixelsrespectively translating the n-pixel algebraically, wherein q ε {1, 2, .. . , n}, and wherein each q-pixel includes (q−1)-faces, (q−2)-faces, .. . , (q-q)-faces; geometrical complexes each being a collection ofq-pixels; q-chains respectively expressing the geometrical complexes inalgebraic form, each q-chain being a linear combination of all theq-pixels of the geometrical complex; in the geometrical complexes,q-cochains which are relations associating quantities related to imagefeatures to the q-pixels and/or faces of said q-pixels; and a coboundarydefining a relation between q-cochains; computing a q-cochain T of afirst of said geometrical complexes as the location of unknowntemperatures; computing a q-cochain H of the first geometrical complexas a global temperature variation; finding a q-cochain ε of a secondgeometrical complex as a global energy variation, as a function of theq-cochain H through a linear transformation; finding the q-cochain ε asa function of the q-cochain T; defining a q-cochain G of the firstgeometrical complex from the q-cochain T through a first coboundaryoperation, transforming the q-cochain G into a q-cochain Q of the secondgeometrical complex, and defining, from the q-cochain Q and through asecond coboundary operation, a q-cochain D of the second geometricalcomplex as a global diffusion; defining a q-cochain S of the secondgeometrical complex as a global source; and establishing a relationbetween the q-cochains ε, D and S.
 45. A computational framework fortwo-dimensional active contour model, comprising: producing an imagesupport including a structure of n-pixels, the image support comprising:q-pixels respectively translating the n-pixel algebraically, wherein q ε{1, 2, . . . , n}, and wherein each q-pixel includes (q−1)-faces,(q−2)-faces, . . . , (q-q)-faces; geometrical complexes each being acollection of q-pixels; q-chains respectively expressing the geometricalcomplexes in algebraic form, each q-chain being a linear combination ofall the q-pixels of the geometrical complex; in the geometricalcomplexes, q-cochains which are relations associating quantities relatedto image features to the q-pixels and/or faces of said q-pixels; and acoboundary defining a relation between q-cochains; computing adisplacement q-cochain D of a first of said geometrical complexes;computing a strain q-cochain S of a second of said geometricalcomplexes, comprising: defining an approximate strain function {tildeover (ε)}(x) as a function of the q-cochain D; expressing the q-cochainS as a function of the approximate strain function and relativepositions of the first and second geometrical complexes; and computing aforce q-cochain F of the second geometrical complex as a coboundary ofthe strain q-cochain S.